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Preface 


The  past  five  years  have  seen  advances  in  the  field  of  computa- 
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in  computer  hardware,  make  future  prospects  very  bright  for  this  field. 
I have  enjoyed  participating  in  the  field  during  this  period  and  would 
recommend  it  to  anyone  who  is  at  all  interested  in  numerics  and  the 
physics  of  fluid  flow. 
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for  his  encouragement,  support,  and  tolerance  of  the  many  alarums  and 
excursions  which  occurred.  I also  wish  to  express  my  appreciation  to 
Dr.  Dale  Anderson  of  Iowa  State  University  for  many  hours  of  fruitful 
discussion  and  for  confirming  that  the  imbedded  shock  instability  en- 
countered was  Indeed  present  in  other  than  the  results  of  my  own  com- 
puter code.  Furthermore,  my  gratitude  goes  to  my  colleagues,  Lt.  Doyle 
Knight  and  Dr.  Joe  Shang,  for  sharing  their  knowledge  of  numerics  and 
fluid  flow  with  me. 

The  credit  for  manuscript  preparation  (and  much  suffering  over 
translating  my  handwriting)  goes  to  Mrs.  Virginia  Cadwell.  Mrs.  Rachel 
Avant  completed  it  after  Virginia  took  a much  deserved  promotion. 

Finally,  thanks  go  to  my  wife,  Linda,  and  children,  Alex  and 
Marcus,  for  tolerating  the  many  absences  over  the  years  (contrary  to 
popular  belief,  computers  only  really  operate  at  night  and  on  weekends) 


The  research  was  performed  between  15  October  1971  and  15  July 
1976  at  the  Theoretical  Aerodynamics  Research  Laboratory,  Aerospace 
Research  Laboratories  (RIP)  and  at  the  Aerodynamics  and  Airframe  Branch 
(AFFDL/FXM),  Wright -Patter son  AFB,  OH. 
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Abstract 


Solutions  were  obtained  for  hypersonic  flow  over  sharp  cones  at 
high  angle  of  attack  which  include  the  viscous  effects  present  in  experi- 
ment. These  solutions  were  obtained  by  integrating  the  Navier-Stokes 
equations  subject  to  a conical  symmetry  assumption.  The  Integration 
technique  used  was  MacCormack's  method  with  boundary  conditions  chosen  to 
match  available  experiments.  The  experimental  conditions  were: 

(1)  Those  of  Tracy  for  a 10°  half  angle  cone  at  angles  of  attack  of  0°, 
8°,  12°,  20°,  and  24°  in  M - 7.95  flow  at  Re  ■ 4.2  x 10**.  (2)  Those  of 
Stetson  for  a 5.6®  half  angle  cone  at  10®  angle  of  attack  in  M - 14.2 
flow  at  Re  • 7.9  x 10^.  (3)  Those  of  McElderry  for  a 6®  half  angle  cone 

at  12®  angle  of  attack  in  M - 6.05  flow  at  Re  • 1.5  x 10^.  A physically 
based  technique  (normal  stress  damping)  was  demonstrated  for  controlling 
starting  transients  and  for  reducing  or  eliminating  numerical  oscilla- 
tions occurring  at  shock  discontinuities  during  the  Integration. 

The  general  features  which  appeared  in  experiment  were  shown  to 
appear  in  the  results  of  the  integration,  including  the  proper  behavior, 
in  laminar  flow,  of  the  viscous  layer  and  the  vortical  singularity.  The 
results  agreed  quite  well  with  the  experimental  data  of  Tracy  (hypersonic 
similarity  parameter  7 =1)  except  for  a discrepancy  in  pitot  pressure 
in  the  viscous  layer  evident  in  a small  region  near  the  leeward  center- 
line  of  the  cone.  The  agreement  with  the  experimental  data  (7  = 3)  of 
Stetson  was  less  adequate.  Surface  pressure  agreement  in  this  instance 
was_  quite  reasonable.  However,  a somewhat  thin  lee  side  viscous  layer 
resulted  in  a calculated  bow  shock  wave  position  27%  closer  to  the  cone 
surface  at  the  lee  centerline.  It  was  concluded  that  the  lee  side 
viscous  layer  discrepancies  at  both  experimental  conditions  were 
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primarily  due  to  lack  of  any  mechanism  In  the  present  technique  to  model 
the  non-conical  flow  near  the  nose  of  the  cone.  A solution  obtained  just 
upstream  of  the  beginning  of  boundary  layer  transition  at  the  experiment- 
al conditions  of  McElderry  agreed  well  with  experiment  when  conically 
projected  far  into  the  turbulent  regime.  The  adequacy  of  the  conical 
symmetry  assumption  is  therefore  Indicated  for  the  turbulent  regime  on 
conical  bodies. 

In  summary,  the  results  show  good  agreement  with  experiment  for 
values  of  the  hypersonic  similarity  parameter  (x  <_  1.0)  and  less 
adequate  agreement  at  higher  x*  Angle  of  attack  limitations  encountered 
in  previous  inviscld  cone  flow  calculations  were  not  encountered  in  the 
present  study. 
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I.  Introduction 


Supersonic  flow  about  cones  has  been  a subject  of  Interest  to 
aerodynamicists  for  many  years.  Cones  comprised  one  of  the  earliest 
attempts  to  streamline  such  items  as  fuel  tanks,  weaponry,  and  fuse- 
lages, with  locally  conical  nose  shapes  still  in  use  on  todayfe  modern 
high-speed  aircraft.  The  missile  era  has  also  seen  the  use  of  cones 
as  a primary  or  partial  shape  for  re-entry  vehicles  in  current  ICBMs. 
This  study  is  directed  toward  the  latter  application  and  is  in  response 
to  the  need  for  more  accurate  aerodynamic  prediction  techniques  in 
this  area. 

The  two  primary  elements  which  dictate  the  need  for  Increased 
prediction  accuracy  are  the  current  high  cost  of  wind  tunnel  testing 
and  the  difficulty  of  simulating  the  re-entry  vehicle  flight  regime. 

The  advent  of  the  high-speed,  large -capacity  computer  in  recent  years 
has  opened  up  a means  to  reduce  cost  both  through  reducing  unit  compu- 
tational cost  and  through  reducing  the  amount  of  experimental  verifi- 
cation required  during  the  systems  acquisition  process.  The  remaining 
task  is  the  research  necessary  to  develop  the  required  numerical  tech- 
niques and  the  verification  of  these  techniques  through  comparison  of 
numerical  results  with  those  obtained  experimentally. 

Particular  solutions  to  the  complete  governing  equations  of  fluid 
flow  are  available  for  very  few  sets  of  boundary  and  initial  conditions. 
It  has  therefore  been  necessary  in  nearly  all  cases  to  reduce  the 
complexity  and  scope  of  the  equations  through  neglect  of  the  viscous 
effects  and/or  by  reduction  of  the  number  of  dimensions  to  be  consid- 
ered in  the  problem.  Prior  to  1973,  most  attempts  to  solve  the  flow 
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about  conical  bodies  used  invlscid  forms  of  the  governing  equations. 

The  two  primary  difficulties  that  have  been  encountered  with  these  solu- 
tions are  the  failure  of  the  inviscid  equations  to  model  properly  the  lee 
side  flow  in  which  viscous  effects  predominate  and  the  occurrence  of 
numerical  instabilities  which  are  encountered  at  high  angle  of  attack. 

If  the  goal  of  reducing  the  amount  of  experimental  verification  is  to  be 
attained,  then  improved  numerical  techniques  must  be  developed  to  over- 
come these  difficulties.  The  present  study  will  present  and  compare  with 
experiment  a technique  which  resolved  both  of  these  difficulties. 

This  report  contains  first  a brief  description  of  the  developments 
in  the  solution  of  flow  over  conical  bodies.  The  solution  approach  used 
in  the  present  ttudy  is  then  presented  in  four  parts.  First  a discussion 
of  the  general  features  of  conical  flows  is  given,  with  particular 
emphasis  on  the  circular  sharp  cone  at  angle  of  attack.  The  governing 
equations  are  then  presented  and  changes  due  to  the  conical  symmetry 
assumption  are  discussed.  The  application  of  MacCormack's  numerical 
integration  scheme  to  the  equations  is  given  in  the  third  j-.art.  Fart 
four  describes  the  details  of  the  finite  difference  mesh  and  the  manner 
in  which  the  integration  was  accomplished.  The  results  of  the  numerical 
integration  are  then  compared  with  experiment  for  selected  conditions 
for  which  previously  published  experimental  data  were  available. 
Conclusions  drawn  from  these  comparisons  are  presented  to  complete  the 
study. 
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II . Background 


The  history  of  the  solution  of  supersonic  cone  flows  has  been 

greatly  affected  by  and  dependent  on  the  development  of  high-speed 

digital  computers.  It  is  generally  agreed  that  the  Navier-Stokes 

equations  in  their  complete  form  provide  sufficient  knowledge  of  the 

physics  of  the  flow  for  non-reacting  gases.  However,  the  scope  and 

conplexlty  of  these  partial  differential  equations  required  great 

simplification  of  the  set  in  order  to  obtain  solutions  prior  to  the 

development  of  the  digital  computer.  These  simplifications  led  to 

poor  agreement  of  the  solution  with  experiment  except  for  restricted 

• • 

conditions.  Although  the  restrictions  have  been  considerably  relaxed 
in  more  recent  studies,  a fully  three  dimensional  solution  of  the 
complete  governing  equations  is  still  not  practical  for  general  config- 
urations. This  section  will  provide  a brief  description  of  the  more 
important  supersonic  cone  flow  solutions  published  prior  to  the 
present  study. 

The  history  of  solution  techniques  for  supersonic  cone  flow  began 
in  1933  with  the  publication  by  Taylor  and  Maccoll  (Ref  1)  of  an 
ordinary  differential  equation  solution  for  axisymmetric  lnviscld  flow 
about  a sharp  cone.  The  next  significant  event  occurred  in  1947  with 
publication  of  tabulated  Taylor-Maccoll  solutions  by  Kopal  and  Asso- 
ciates (Ref  2) . These  tabulations  were  extended  to  include  angle  of 
attack  in  Refs  3 and  4.  Ferrl  (Ref  ’>)  obtained  first  order  solutions 
for  cones  at  small  angle  of  attack  in  which  he  noted  the  properties  of 
the  vortical  layer  about  the  cone  and  showed  the  existence  of  a vortical 
singularity  in  the  lee  side  flow.  With  this  publication  by  Ferri,  all 
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important  features  (except  for  the  lee  side  Imbedded  shock  waves)  of 
the  invlscid  conical  flow  fields  were  known  and  defined.  Since  Ferri's 
solution  was  demonstrated  only  for  ot/@c  £ 1.6,  it  remained  then  to 
Increase  the  angle  of  attack  for  which  a solution  could  be  obtained  and 
to  account  for  imbedded  shock  waves  and  viscous  effects  present  in 
experiment. 

Improvement  of  the  digital  computer  in  the  middle  1960's  allowed 
the  solution  of  the  non-linear  Euler  equations  in  more  complete  form. 
Babenko,  et  al  (Ref  6)  obtained  smooth-body  solutions  for  the  lnviscld 
Euler  equations  in  supersonic  flow  without  further  simplification. 

This  work  is  an  early  example  of  the  use  of  finite  difference  tech- 
niques to  integrate  the  three-dimensional  flow  equations  in  their  non- 
linear form.  Another  example  of  this  approach  is  the  work  of 
Bohachevsky  and  Rubin  (Ref  7) . 

The  finite  difference  scheme  used  in  the  present  work  can  be 
traced  back  to  that  of  Lax  and  Wendroff  (Ref  8).  This  second  order 
scheme,  when  applied  to  the  flow  equations  case  in  conservation  law 
form  allows  shock  discontinuities  to  form  without  special  provisions 
during  the  integration  process.  MacCormack's  (Ref  9)  variant  of  the 
Lax-Wendroff  scheme  as  applied  by  Kutler  (Ref  10)  extended  the  angle  of 
attack  range  for  which  invlscid  sharp  cone  flows  could  be  obtained. 
Further  applications  (Refs  11-13)  of  this  technique  to  invlscid  flows 
were  made  by  Kutler,  et  al  to  allow  for  real  gas  effects,  varying 
geometries,  and  passing  through  shock  fronts.  However,  the  inaccuracies 
due  to  neglect  of  viscous  effects  were  not  effectively  dealt  with  until 
the  work  of  Lin  and  Rubin  (Refs  14,13)  and  Lubard  and  Helliwell  (Ref  16). 
Lin  and  Rubin  solved  a boundary  layer  equation  set  modified  to  account 
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for  centrifugal  force  and  cross  flow  diffusion  in  the  weak  interaction 
region  and  a parabolic  set  in  the  hypersonic  tip  region  to  obtain 
solutions  including  viscous  effects  up  to  twice  the  cone  half  angle. 
However,  the  former  approach  requires  input  of  surface  pressure  from 
experiment  or  other  source.  Lubard  and  Helliwell  utilized  an  Implicit 
finite  difference  technique  similar  to  that  of  Lin  and  Rubin  to  solve 
the  flow  equations  subject  to  a parabolic  assumption  for  the  stream- 
wise  shear  stress  terms.  This  space  marching  technique  required  an 
accurate  initial  surface  for  the  integration  to  proceed.  The  lack  of 
such  a surface  required  that  the  integration  be  started  at  zero  angle 
of  attack  and  that  the  angle  of  attack  be  gradually  increased  until 
the  desired  value  is  reached.  The  solution  must  then  be  carried  far 
enough  downstream  for  relaxation  of  the  effects  of  this  procedure  to 
occur.  The  comparison  with  experiment  for  this  technique  at  a given 
physical  point  on  the  cone  surface  is  therefore  uncertain. 

As  noted  in  this  brief  survey,  solutions  of  viscous  flow  over 
cones  have  only  recently  been  attempted  (1973).  The  present  study 
seeks  to  illustrate  a technique  which  Includes  viscous  effects  in  cone 
flow  calculations  and  removes  limitations  present  in  the  inviscid 
calculations.  The  details  of  this  technique  are  given  in  the 
following  section. 
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III.  Approach 


This  section  presents  the  details  of  the  solution  approach  used  in 
the  present  study.  A brief  description  of  the  general  features  of  coni- 
cal flow  is  presented  in  part  one.  This  is  followed  by  a discussion  in 
part  two  of  the  governing  equations  of  fluid  flow  and  the  changes  in 
these  equations  which  result  from  the  conical  symmetry  assumption.  Part 
three  describes  the  application  of  MacCormack's  numerical  integration 
scheme  to  the  equation  set  obtained  in  part  two.  Part  four  then  gives 
the  details  of  the  solution  procedure. 

Description  of  Conical  Flow 

As  noted  in  the  introduction  to  the  present  study,  supersonic  flow 
about  conical  bodies  has  been  of  considerable  interest  to  aerodynamlcis ts . 
Supersonic  inviscld  conical  flows  appear  at  first  inspection  to  be  fully 
three-dimensional.  However,  they  are  in  reality  only  two  dimensional 
with  no  gradients  in  any  quantity  occurring  along  the  third  dimension. 

A more  detailed  description  of  these  flows  is  presented  here  to  promote 
understanding  of  their  unique  properties. 

Truly  conical  supersonic  flows  are  steady  state  solutions  to  the 
Euler  equations  (Navier-Stokes  equations  with  viscous  terms  neglected) 
for  a particular  set  of  boundary  conditions.  These  conditions  are  a 
uniform  supersonic  outer  flow  and  a body  generated  by  rays  passing 
through  a common  vertex.  To  ensure  that  the  flow  is  completely  conical, 
the  additional  condition  of  a bow  shock  wave  attached  to  the  vertex  of 
the  body  must  be  imposed.  Examples  of  bodies  which  generate  these  flows 
are:  (1)  circular  cones,  (2)  elliptical  cones,  (3)  delta  wings, 

(4)  delta  wlng/conical  body  combinations,  and  (5)  axial  corners. 
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Many  variations/combinations  of  these  bodies  are  possible  which  produce 
conical  flows.  Also,  conical  bodies  which  appear  as  nose  shapes  in  more 
general  configurations  will  produce  locally  conical  flows.  The 
characteristics  of  these  flows  are  examined  below. 

The  dominant  feature  of  conical  flow  is  that  all  fluid  quantities 
are  constant  along  rays  passing  through  the  vertex  of  the  conical  body. 

By  noting  that  these  rays  are  equivalent  to  the  spherical  radius  r with 
origin  at  the  cone  vertex,  conical  flow  can  be  mathematically  described 
by  stating  that  derivatives  of  all  fluid  quantities  with  respect  to  r are 
Identically  zero.  The  physical  Implications  of  this  statement  are 
illustrated  in  Fig.  1.  This  figure  shows  a plane  cut  of  an  axlsymmetrlc 
supersonic  cone  flow.  A conical  flow  exists  when  the  flow  quantities  at 
point  A (along  streamline  are  identical  to  those  at  point  B (along 
streamline  ip  ).  The  flow  at  both  point  A and  point  B can  be  completely 

D 

described  in  terms  of  the  single  angle  8 for  this  two  dimensional  cut. 
What  is  apparently  a two  dimensional  flow  is  then  in  reality  only  one 
dimensional  (sometimes  referred  to  as  1-1/2  dimensional,  since  two 
velocity  components  are  still  present) . 

Although  somewhat  more  difficult  to  visualize,  conical  bodies  at 
angle  of  attack  in  supersonic  flow  also  produce  conical  flow  fields.  By 
analogy  with  the  above  discussion,  these  flows  can  also  be  described 
completely  in  terms  of  the  independent  variable  6 and  a circumferential 
angle  <j>  (these  angles  are  illustrated  in  Fig.  2).  All  natural  features 
of  the  flow;  such  as  imbedded  shocks,  bow  shocks,  slip  lines,  etc.;  form 
surfaces  composed  of  rays  passing  through  the  conical  vertex.  An  example 
sketch  of  this  flow  is  shown  in  Fig.  3.  The  streamlines  shown  in  Fig.  3 
show  the  manner  in  which  the  flow  crosses  the  bow  shock,  expands  around 
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the  body,  and  then  is  turned  back  along  the  cone  by  the  imbedded  shock. 
It  is  instructive  to  view  flow  features  on  a spherical  surface  through 
the  field. 

Figure  4 illustrates  the  approximate  cross  flow  streamline  pattern 
on  the  spherical  surface  and  shows  relative  location  of  the  shock 
waves.  The  point  at  which  the  cross  flow  streamlines  converge  is  the 
vortical  singularity  first  noted  by  Ferri  (Ref  5) . The  streamlines 
all  have  different  values  of  entropy  which  Implies  that  density  and 
entropy  are  discontinuous  at  the  singularity  with  pressure  being 
continuous.  This  singularity  rests  on  the  cone  surface  (in  fact, 

Melnik  (Ref  17)  has  shown  the  possibility  of  two  singularities  on  the 
surface)  at  low  angles  of  attack  and  lifts  off  as  angle  of  attack 
increases. 

In  summary,  conical  flows  can  be  completely  described  in  a 
spherical  coordinate  system  by  only  two  independent  coordinates  (0,  <J>) . 
All  natural  features  of  the  flow  are  also  conical  surfaces.  These 
facts  tend  to  render  a very  complex  flow  more  amenable  to  solution  by 
presently  available  techniques. 

The  Fluid  Flow  Equations 

In  this  section,  the  basic  equations  of  fluid  flow  are  presented 
in  complete  form  and  the  conical  symmetry  assumption  is  applied.  The 
resulting  equations  allow  integration  to  take  place  on  a single 
spherical  surface.  The  non-dimenslonalization  of  the  equations  and 
the  physical  meaning  of  the  assumptions  are  discussed. 
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The  Navier-Stokes  equations  which  describe  £low  of  a perfect  gas, 
can  be  written  in  conservation  law  form  (Ref  18)  for  a spherical 
coordinate  system  as  follows: 


3D  3E  3F  3G  + 

3t  + 3r  + 30  + 3<J>  + K " 0 


(1) 


where 


rz  sin6 


1 • r2  sin0 


P 

pu 
pv 
pw 
pe 

pu 

pu2  - an 
puv  - T l 2 
puw  - T 1 3 


pue  + qr  — U0; i 


- VT  j 2 - WT 1 3 


r sind 


pv 

pUV  - T 2 1 
PV2  - a22 
PWV  - T 2 3 


pve  + q0  - ut2  i 


- V02  2 - WT  2 3 
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puw  - T31 


pvw  - 1 3 2 


PW  - C7  3 3 


pwe  + - ui3i 


- VT  32  - WO  3 3 


H ™ r sin0 


-(pV2-022)  - (PW2  -O33) 
-cot0  (pw2— cj  33)  + ( pUV— T 12) 
COt0  ( PVW-T  23)  + (pUW-Tis) 


The  stress  tensor  terms  and  will  be  defined  in  detail  later. 

Solution  of  the  equation  set  in  the  form  shown  above  would  require 
very  large  amounts  of  machine  time  even  for  simple  aerodynamic  config- 
urations. This  has  led  researchers  in  the  past  to  simplify  the 
equation  set  through  neglect  of  the  stress  terms  (inviscld  flow) 
and/or  by  reducing  the  number  of  dimensions  to  be  considered  in  the 
problem.  One  class  of  flows  which  has  been  extensively  examined, 

(Refs  10,19,20)  through  neglect  of  the  viscous  terms,  is  that  of 
inviscld  conical  flow. 

A conical  flow,  as  noted  in  the  previous  section,  can  be 
described  as  an  inviscld  steady  flow  in  which  all  flow  quantities  are 
constant  along  rays  passing  through  the  vertex  of  the  conical  body. 

If  a body  fixed  spherical  coordinate  system  (Fig.  2)  centered  at  the 
vertex  of  the  conical  body  is  used  to  describe  the  flow,  then  all 
spherical  surfaces  must  have  the  same  vector  and  scalar  values  of  the 
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flow  quantities  for  a given  (6,  <j>)  point  on  each  surface.  Therefore 
all  derivatives  of  flow  quantities  with  respect  to  the  spherical 
radius  (r)  of  these  spherical  surfaces  from  the  origin  must  be  zero. 
This  has  the  effect  of  reducing  the  number  of  independent  variables  in 
the  problem  by  one. 

Examination  of  experimental  studies  of  supersonic  flow  over 
conical  bodies  (Refs  21-23)  reveals  that  these  flows  exhibit  approxi- 
mate conical  behavior  downstream  from  the  nose  region  even  though 
relatively  large  viscous  regions  exist.  Cross  (Ref  21)  determined 
that  the  viscous  layer  growth  on  the  lee  side  of  a delta  wing  in  super- 
sonic flow  was  essentially  conical.  The  oil  flow  separation  traces  for 
sharp  cones  in  the  experimental  study  of  Stetson  (Ref  22)  are  approxi- 
mately straight  (but  the  conical  vertex  of  these  traces  is  displaced 
downstream  by  nose  effects).  Therefore,  in  concert  with  an  idea  first 

broached  by  Anderson  (Ref  24)  for  axial  corner  flow,  the  conically 

g 

symmetric  flow  approximation  (-^  = 0)  is  applied  to  all  terms  in 
Eq  (1). 

The  resulting  equation  set  is: 


5D  3F 
3t  30  + 3$ 


+ H" 


0 


(2) 


where : 

(a)  D,  F,  and  £ are  unchanged  except  in  the  definition  of  the 
stress  terms. 
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3 

[ 1 

3r 

r2sin0 

0 

au 

t2i 

T31 

uan  + vt  1 2 + wr  1 3 


(c)  t'  - H + ~ + J r2sin0 


or  H""  ■ r sin©  | 2 pu 

-(pv2-022)  - (pV2-C3j)  + 2pu2-an 
-COt0  (pW2  — C 33)  + 3puv-2ti2 
cofc0  (pvw—Tj  3)  + 3puw-2ti3 
2pue  - uan  “ vT  12  - wr  1 3 
These  equations  are  non-diraensionalized  as  follows,  with  dimensional 
quantities  denoted  by  primes: 


V 

max 


V 

max 


V 

max 


P 


P 


e “ e 
V 

max 


t - t 

r7v 


T - T' 


max 


Re  ■ Vmax  Pt<»r 


V 

max 


lc 


-Hi 

2y 


(3) 


Equation  3 indicates  a notational  change  from  this  point  onward. 
Hereafter,  units  will  be  specified  where  dimensional  quantities  are 
used.  The  quantities  used  for  non-dimens  ionalizatlon  are  due  to 
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Kutler  (Ref  10) . It  is  significant  that  when  time  is  non-dimen3ion- 

alized  by  the  parameter  r/V  , the  r dependence  in  Eq  (2)  is  contained 

max 

in  the  Reynolds  number.  The  net  effect  is  that  the  calculation  is 
carried  out  at  a single  (8,  <f>)  spherical  surface  with  the  distance  of 
this  surface  from  the  cone  apex  determined  by  the  Reynolds  number. 
Therefore,  all  spherical  radius  scaling  is  now  contained  in  the 
Reynolds  number  alone. 

The  non-dimensional  stress  terms  for  this  equation  set  are 
defined  as  follows: 


t12 

Tis 

t23 

a 


n 


On 

<*22 


T!1  ' IS  <!f  - v> 

Tl‘  " RTTISe'l?  ’ ” 8ln0) 

T”  ' STSIHe  [I?  <w  sl"9)  -2  w c°s9  + If] 

Ve  [2  " s1”9  +!e  (v  ain9)  + If  j 


2 

3Re 


- kp  + 


sin0 


, . 2 ,3v  . . , n 

' k»  + iu(96  + u)  + il «S 


- kp  + 


r9w 


Re  sln6  9<j> 


(t-t  + u sin  + v cos0)  + 


sln6 


(4) 


<733 


The  heat  transfer  terms  are  defined  as  follows: 


‘ 1 3T  • 1 3T  ... 

q6  " “ 2RePr  36  » q<|>  " “ 2RePr  3<J>  W 

VC 

where  Pr  ■ — ^ 

At  this  point,  it  is  prudent  to  examine  further  the  physical 
Implications  of  the  approximations  applied  to  the  above  equation  set. 

It  should  first  be  noted  that  the  only  approximations  now  Inherent  to 
the  equations  are  that  a perfect  gas  is  required  and  that  conically 
symmetric  flow  is  assumed  at  all  points  in  the  field.  Figure  5 
depicts  a cross  section  through  an  axlsymmetric  cone  flow  with  a some- 
what exaggerated  boundary  layer  thickness  shown  by  the  dashed  line. 

The  thesis  is  that  at  a given  calculation  surface,  the  viscous  layer 
thickness  will  be  properly  scaled  by  the  Reynolds  number  based  on  the 
spherical  radius  to  that  surface.  The  only  error  then  would  be  in  the 
local  gradients  of  the  flow  quantities  and  would  be  of  the  order  of 

the  angle  between  the  spherical  radius  and  the  edge  of  the  boundary 

layer  at  the  calculation  surface.  The  same  statement  could  be  made 
concerning  a calculation  at  any  point  along  the  cone  downstream  from 
the  non-conical  nose  region.  Therefore,  since  the  growth  of  the 
viscous  layer  thickness  will  not  in  general  be  linear,  the  assumption 
in  reality  only  Implies  locally  conical  flow. 

In  order  to  test  the  above  hypothesis,  calculations  were  carried 
out  using  the  above  equation  set  and  the  method  set  forth  in  the 
following  sections  applied  to  axlsymmetric  cone  flow.  The  cone  surface 

pressure  resulting  from  these  calculations  is  shown  in  Fig.  6 and 

compared  with  hypersonic  weak  interaction  theory  for  axlsymmetric  cone 
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flow  (Ref  29) . The  agreement  shown  in  Fig  6 between  the  calculations 
and  the  weak  Interaction  theory  Is  quite  good  with  divergence  noted 
(as  expected)  toward  the  nose  of  the  cone.  In  order  to  fully  verify 
the  hypothesis,  two  dimensional  calculations  compared  with  appropriate 
experimental  data  are  required.  The  remainder  of  this  report  will 
present  the  results  of  calculations  carried  out  at  selected  conditions 
for  which  experimental  data  are  available.  Comparisons  are  made  and 
conclusions  are  drawn  concerning  the  adequacy  of  the  hypothesis. 

The  Numerical  Integration  Scheme 

This  section  discusses  the  choice  of  MacCormack's  finite  differ- 
ence scheme  as  the  integration  method  for  this  problem.  The  scheme  is 
illustratively  applied  to  the  one-dimensional  wave  equation  as  an 
example.  The  equations  to  be  integrated  are  cast  in  MacCormack's 
predictor-corrector  form,  and  the  finite  difftrence  representation  of 
the  derivatives  is  shown. 

The  utility  of  MacCormack's  (Ref  9)  finite  difference  scheme  for 
the  solution  of  inviscld  supersonic  flows  has  been  demonstrated  in 
numerous  studies  (Refs  10-13)  . More  recent  studies  have  also  shown  it 
to  give  excellent  results  for  two-dimensional  viscous  flows  with  sepa- 
ration (Refs  25-27) . The  scheme  can  be  described  as  a second  order 
accurate  predictor-corrector  sequence  for  the  integration  of  partial 
differential  equations.  MacCormack's  scheme  is  a variant  of  the  Lax- 
Wendroff  scheme  (Ref  10)  and  can  be  shown  (as  applied  to  the  linear 
wave  equation)  to  reduce  to  Lax-Wendroff  form  when  the  predictor  is 
substituted  into  the  corrector.  Since  second  order  difference  schemes 
generally  give  better  results  with  fewer  mesh  points  than  first  order 
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schemes,  MacCormack's  proven  differencing  scheme  Is  chosen  for  the 
present  study. 

If  an  Euler  predictor  followed  by  a modified  Euler  corrector  is 

8U  3U 

applied  to  the  1-D  wave  equation  + “ 0; 


Predictor: 

*1+1 

■ U1  - <Ux>i 

(6) 

Corrector: 

Ui+1 

- j[Di + K+i  - c4t 

<»*>i+i] 

(7) 

MacCormack's  scheme  results  when  forward  first  spatial  differences  are 
used  in  the  predictor  and  backward  first  differences  in  the  corrector. 
The  predictor  corrector  sequence  is  then  (with  i representing  time 
location  and  j spatial  location) : 


Predictor:  lf*+1  - uj  - ^ (D*+1  - U^)  (8) 

Corrector:  uf1  - ± g (9J+1  - »£)]  (9) 


The  two  equations  collapse,  when  the  appropriate  substitutions  are 
made,  to  the  following  Lax-Wendroff  form: 


Ui+1  - U*  - - t - (U1  - U1  ) + 

J 3 - Jt  V j+1  J-r 


:2At2 

"te7 


- 2U1 + »J.i> 


(10) 


The  leading  error  terms  which  are  truncated  through  use  of  this  finite 
difference  equation  representation  are  of  the  form  (with  the  derivatives 
evaluated  at  an  undetermined  point  within  the  range  of  interest) 


cAtAx2  cAt2Ax2  At3 

+ ~ Ux3  ' V f T~  Ut3 
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which  Illustrates  that  when  At  and  Ax  are  of  the  same  order,  the 
estimated  truncation  error  must  be  of  0(AtJ)  and  the  method  can  be 
said  to  be  second  order  accurate.  The  stability  of  MacCormack's 
scheme  Is  examined  in  Appendix  C. 

When  the  scheme  is  applied  to  Eq  (2),  the  resulting  MacCormack's 
predictor-corrector  steps  are  as  follows  (the  used  previously  to 


indicate  vector  quantities  is  omitted  here  for  clarity); 
Predictor  - F^)  - 


Corrector 


- G,1.  ) - At  E'  . 

A4>  j.k+1  j,k'  j,k 


(id 


(12) 


where  ^ indicates  that  the  flow  quantities  are  evaluated  at  the 

predictor  level  and  D*  , Implies  D(lAt,  jA0,  kA<J>) . The  presence  of 

J 

the  H"  matrix  in  these  equations  indicates  that  they  are  in  so  called 
"weak"  conservation  form  (Ref  28) . Numerical  results  in  the  present 
study  reveal  that  no  significant  error  in  total  temperature  occurs 
across  shock  transitions  through  use  of  this  form  of  the  equations 
(implying  that  energy  must  have  been  properly  conserved). 

Since  Eq  (2)  contains  derivatives  in  the  stress  terms  which 
remain  as  an  Integral  part  of  the  matrix  terms,  a value  of  the  deriva- 
tives must  be  obtained  at  each  of  the  mesh  points  used  in  the  primary 
differencing.  These  values  must  be  obtained  in  a manner  which  vrill 
maintain  both  consistency  and  accuracy  of  the  differencing  equation. 
This  was  done  in  the  present  study  by  evaluating  the  interior  derivative 
terms  with  respect  to  the  same  Independent  variable  using  first  order 
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difference  approximations  of  opposite  sense  to  the  exterior  differences. 
An  illustration  follows  for  the  predictor. 

Consider  the  vector  F which  consists  of  functions  of  the  flow 
quantities  and  their  derivatives.  The  finite  difference  representa- 
tion of  -Ijr  is 


m * i_  a _ t ) 

36  A0  'Fj+l,k  Fj,k; 


(13) 


. 3u 

Assume  that  the  vector  F contains  the  derivative  values  of  which 

are  therefore  required  at  mesh  points  j,k  and  j+l,k  in  the  difference 

3u  1 

representation.  These  values  are  obtained  by  = (u..,  , - u ) 

dw  Zav  j+ljR  j 

at  J+l,k  and  (u^  - u^  at  j,k.  The  resulting  finite 

32u 

difference  representation  for  is  the  standard  second  order  accurate 
difference  centered  at  j . 


3?. 


When  a cross  derivative  term  is  involved  (say  a ^ term  contained 


3F 


in  -gg-)  , second  order  centered  first  differences  are  used.  With  gg 

3u  1 

represented  as  above,  ^ ^ (uj+l,k+l  “ “j+l.k-15  at  J+1’k  *** 

3u  ^ 

(u.  , - u.  . ,)  at  1,k.  This  results  in  the  standard  second 

o<p  A<{>  j,k+l  J ,k-l 

32u 

order  accurate  second  difference  for  the  difference  centered 


3<J>30 


at  j+*s. 


This  procedure  is  also  used  in  the  corrector  step  with  all  differ- 
encing reversed  in  the  6 direction.  The  resulting  differences,  when 
combined  with  the  predictor  into  the  Lax-Wendroff  form,  are  correct 
and  centered  at  j . 

This  essentially  completes  the  application  of  MacCormack's  scheme 
to  Eq  (2).  It  should  be  apparent  that  throughout  the  analysis  of 
MacCormack's  scheme,  careful  prescription  of  the  first  order  differen- 
ces la  the  individual  predictor  and  corrector  results  in  a second  order 
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accurate  scheme  for  the  entire  cycle.  The  details  of  the  integration 
procedure  are  described  in  the  following  section. 

The  Solution  Procedure 

In  this  section,  the  details  of  the  integration  procedure  are 
discussed.  The  finite  difference  mesh  is  first  described,  followed  by 
the  initialization  used  for  the  initial  value  problem  and  a brief 
description  and  discussion  of  the  physical  boundary  conditions. 
Finally,  the  convergence  criteria  for  determining  when  the  sought-for 
steady  state  solution  has  been  reached  is  described. 

As  stated  earlier  in  this  section  the  calculation  takes  place  on 
a spherical  surface  located  at  a distance  r from  the  vertex  of  the 
cone.  The  distance  r is  determined  by  the  distance  along  the  cone  at 
which  correlation  with  experiment  is  desired.  It  then  appears  as  the 
characteristic  length  in  the  Reynolds  number  and  scales  the  viscous 
effects  at  the  calculation  surface.  The  integration  procedure,  as 
with  all  finite  difference  techniques,  requires  values  of  the  fluid 
flow  quantities  to  be  known  at  discrete  points  labeled  as  mesh  points. 
For  the  present  study,  these  points  are  arranged  on  the  (0,  0)  calcu- 
lation surface  as  shown  in  Fig.  2.  The  mesh  points  are  evenly  spaced 
in  the  0 and  <f>  directions  with  A0  + A$.  It  should  be  noted  that  con- 
siderable advantage  for  viscous  calculations  can  sometimes  be  obtained 
by  varying  mesh  point  spacing  to  cluster  points  near  the  surface. 
However,  only  uniform  size  mesh  is  considered  here  since  the  present 
study  is  primarily  concerned  with  proof  of  an  untried  assumption  for 
the  flow  equations.  Therefore,  the  additional  complexity  involved  in 
obtaining  the  uneven  mesh  spacing  was  not  desirable  and  was  not  made 
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a part  of  this  study.  The  mesh  spacing  used  was  designed  to  give 
entirely  adequate  results  on  the  lee  side  of  the  cone  in  the  large 
viscous  regions  present  there. 

The  finite  difference  mesh  is  Initialized  to  free  stream  values 
of  the  flow  quantities  except  for  surface  points  at  which  the  three 
velocity  components  are  set  identically  to  zero.  This  is  physically 
equivalent  to  bringing  the  cone  from  rest  to  the  free  stream  Mach 
number  Instantaneously.  As  might  be  expected,  the  initial  transients 
associated  with  this  procedure  are  quite  large  and  require  numerical 
damping  at  the  higher  angles  of  attack.  This  is  not,  however,  consid- 
ered detrimental  since  the  damping  is  required  continuously  for  most 
of  the  solutions  obtained  for  reasons  which  will  be  discussed  later. 
Also,  a comparison  between  final  results  obtained  with  the  impulsive 
start  and  with  a solution  obtained  by  changing  slightly  the  boundary 
condition  representation  (thereby  perturbing  the  flow)  of  a converged 
solution  revealed  no  essential  difference  between  the  two  solutions. 

The  boundary  conditions,  which  are  maintained  throughout  the 
calculation,  are: 

(1)  All  velocities  identically  zero  at  the  surface. 

(2)  Surface  temperature  constant  and  equal  to  the  value  for  the 
experimental  case. 

(3)  Free  stream  values  of  all  flow  quantities  are  maintained  at 
the  outer  edge  of  the  mesh.  Since  the  bow  shock  wave  was  very  close 
to  the  cone  surface  on  the  windward  side  for  the  high  a cases,  an 
ellipse  was  used  to  provide  a dummy  outer  edge  of  the  mesh  which  would 
be  close  to  the  bow  shock  and  therefore  save  computer  time. 
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(4)  Lateral  symmetry  0f  the  £low  quantities  Is  maintained  across 
the  $ ■ 0°  and  $ ■ 180°  lines. 

The  values  of  pressure  and  density  at  the  surface  are  determined  during 
the  Integration.  A discussion  of  this  matter,  plus  a more  complete 
description  of  the  boundary  conditions,  appears  in  Appendix  A. 

In  this  study,  the  predictor  and  corrector  are  swept  through  the 
(<(»,  6)  mesh  In  turn,  beginning  at  the  <J>  ■ 0°,  0 * 0c  point.  A complete 
traverse  of  the  predictor  and  corrector  constitutes  one  time  step. 

Step  by  step,  the  procedure  is  as  follows: 

(1)  The  predictor  is  swept  through  the  mesh  from  0 * 0£  to  0 ■ 9® 

and  from  <J>  ■ 0°  to  180°. 

• • 

(2)  The  D matrix  is  decoded  to  obtain  updated  values  of  the  flow 
quantities. 

(3)  The  corrector  is  swept  through  the  mesh,  using  the  updated 
flow  quantities  where  appropriate. 

(4)  The  b matrix  is  decoded  to  give  values  for  the  flow  quanti- 
ties at  the  end  of  the  time  step. 

(5)  Boundary  conditions  are  applied  where  necessary  during  the 
integration. 

As  implied  in  the  stated  outer  boundary  condition,  no  special 
provisions  are  included  to  handle  the  bow  shock  wave  or  any  imbedded 
shock  waves  which  may  arise.  These  shock  waves  are  "captured"  through 
use  of  the  conservation  form  of  the  governing  equations  (Refs  10,28). 
This  procedure  is  considered  adequate  through  use  of  numerical  damping 
for  reducing  spurious  oscillations  around  the  shock  waves  (Appendix  B) . 

A criteria  for  determining  when  the  steady  state  has  been  reached 
(i.e.  "convergence")  is  considered  an  integral  and  necessary  part  of 
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this  study.  Previous  unpublished  experience  by  the  author  with 
lnviscid  flows  reveals  that  instabilities  in  the  calculation  can  arise 
even  after  changes  in  the  flow  field  pressure  become  small  as  seen  on 
a plotter  or  CRT  display  device.  Therefore,  a stringent  convergence 
criteria  is  used  in  the  present  study  in  order  to  insure  that  the 
solution  is  stable  and  that  the  results  obtained  are  repeatable.  The 
convergence  criteria  that  was  used  required  computation  of  differences 
in  pressure,  density,  velocity,  and  energy  between  time  steps  and 
stopped  the  calculation  when  these  differences  reached  the  fifth  signi- 
ficant figure  for  all  points  in  the  flow  field.  At  the  time  step 
sizes  used  for  runs  reported  here,  convergence  of  the  solution  typically 
occurred  at  the  physical  time  required  for  the  outer  flow  to  move  4 to 
5 times  the  length  from  the  nose  to  the  calculation  surface.  By 
reducing  the  number  of  significant  figures  input  to  the  convergence 
routine,  it  was  found  that  the  fifth-place  criteria  gave  run  times 
nearly  twice  as  long  as  would  have  been  required  to  obtain  values  of 
surface  pressure  accurate  to  engineering  standards  (Changes  of  less 
than  3%  occurred  in  surface  pressure  during  the  last  50Z  of  the  run) . 

The  numerical  stability  of  the  present  calculations  is  therefore  assured. 

The  present  study  did  not  make  continuous  use  of  a stability 
limit  criteria  to  determine  allowable  step  size.  Instead,  the  time 
step  size  for  each  computer  run  was  frozen  at  a value  which  would 
Insure  no  difficulty  with  normal  stress  damping  (described  in  Appendix 
B) . Based  on  experience  gained  in  the  study,  a stability  criteria  was 
postulated  and  confirmed  by  numerical  experiment.  This  entailed 
repeat  of  a previous  run  with  time  step  size  automatically  determined  by 
applying  the  criteria  noted  below  at  each  time  step.  This  stability 
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limit  criteria  is  described  in  aetail  in  Appendix  C.  The  criteria  is 


summarized  here: 

At  - m MIN  <Atcfl,  AtnsJ,  Atyl,c>  (14) 


In  this  instance,  MIN  implies  the  minimum  value  of  the  quantity 
found  by  searching  all  mesh  points  in  the  flow  field.  Although  m is 
theoretically  limited  to  1.0  (the  CFL  condition),  it  was  found  that  a 
value  of  m * 1.2  could  be  used  once  the  Initial  strong  transients 
were  past.  Instability  of  the  solution  occurred  at  m ■ 1.4. 

The  operation  sequence  used  in  the  computer  code  was  the 
following : 

(1)  READ  input  data 

(2)  Initialize  mesh 

(3)  DO  time  steps,  1 to  input  number 

(4)  Calculate  local  Reynolds  number 

(5)  Calculate  local  damping  coefficients 

(6)  DO  predictor  for  all  but  free  stream  boundary  mesh  points 

(a)  Calculate  coordinate  metrics 

(b)  CALL  boundary  value  subroutine 

(c)  Define  finite  difference  representations 

(d)  CALL  vector  load  subroutine  (creates  F,  G,  and  H"  in 
Eq  (2)  according  to  the  requirements  of  (c)  above) 
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(e)  CALL  integrations  subroutine  (creates  predicted  value 
of  D at  new  time  level  by  application  of  Eq  (11) 

(f)  CALL  D vector  solver  subroutine  (decodes  flow  quanti- 

•> 

ties  from  calculated  D vector  at  new  time  level) 

(7)  CALL  symmetry  boundary  condition  subroutine 

(8)  DO  corrector  for  all  but  free  stream  boundary  mesh  points 
(This  is  an  exact  duplicate  of  step  (6)  except  for  changes  in  finite 
difference  representations  and  flow  quantity  time  levels.  The  code 
used  the  same  subroutines  for  both  predictor  and  corrector  sequences.) 

(9)  CALL  symmetry  boundary  condition  subroutine 

(10)  Test  for  convergence.  If  not  converged,  return  to  step  (3) 
and  continue  until  convergence  criteria  is  met  or  the  input  number  of 
time  steps  is  completed. 

(11)  Store  results 

(12)  Print  and  plot  results  of  calculation 

(13)  STOP 

This  sequence  of  operations  may  be  used  to  Integrate  (by  use  of 
MacCormack's  scheme)  other  three  coordinate  versions  of  the  flow 
equations  by  appropriate  changes  in  steps  (6)  (b) , and  (d) . The 
restrictions  are  that  weak  conservation  form  be  used  and  that  the 
independent  variable  of  integration  be  time  or  a coordinate  direction 
in  which  M >.  1. 

This  completes  the  description  of  the  approach  used  in  the  present 
study.  The  following  section  describes  and  compares  with  experimental 
data  the  numerical  results  obtained  through  application  of  this 
approach. 


24 


IV.  Numerical  Results 


This  section  presents  a comparison  between  results  obtained  with  the 
present  technique,  the  numerical  technique  of  Lubard  and  Helliwell,  and 
selected  wind  tunnel  experiments.  The  first  comparison  is  with  the 
experimental  study  of  Tracy  (Ref  30).  This  study  has  been  a standard  of 
comparison  for  other  techniques  (Refs  14-15)  and  Is  very  well  documented 
in  the  open  literature.  The  second  comparison  is  with  the  experimental 
study  of  Stetson  (Ref  22)  and  the  third  with  that  of  McElderry  (Ref  31) . 

The  validity  of  any  heretofore  untried  assumption  applied  to  the 
governing  equations  can  only  be  determined  through  comparison  of  the 
results  obtained  by  calculation  with  results  of  experiment.  To  test  the 
validity  of  the  present  technique,  experimental  data  were  chosen  which 
had  been  published  and  which  provided  a Mach  number  range  of  six  to 
fourteen.  The  experimental  cases  chosen  were  for  sharp  cones  within  the 
limits  reasonably  obtainable  by  standard  machining  techniques  (0.001  to 
0.002  in.  nose  radius).  Reynolds  numbers,  except  for  that  of  McElderry, 
were  in  the  range  for  which  laminar  flow  could  be  expected.  Transition 
onset  in  McElderry' s data  was  at  approximately  Re  * 2.4  x 10^.  The 
cases  all  have  flow  field  probe  data  available  for  selected  Reynolds 
numbers  and  angles  of  attack.  The  remainder  of  this  section  will 
present  a detailed  comparison  of  the  results  of  the  present  calculation 
with  the  chosen  experimental  data. 
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Tracy  * s Data 


Tracy  (Ref  29)  conducted  an  experimental  study  using  a 10°  half 
angle  sharp  (D. 002  nose  radius)  cone  in  M - 7.95  flow.  The  wind  tunnel 
working  fluid  was  air  with  total  pressure  and  temperature  of  259.3 
psla  and  1360°R  for  the  runs  of  interest.  The  resulting  free  stream 
Re  was  1.25  x 10^/Ft.  Data  were  taken  at  angles  of  attack  of  0°  to  24° 
in  4°  increments.  Measurements  of  surface  pressure  at  x ■ 4.0  in.  and 
surveys  of  picot  pr*»p»ure  at  x ■ 3.45  in.  (and  perpendicular  to  the 
cone  surface)  were  selected  for  comparison.  The  configuration  of  the 
model  and  instrumentation  prevented  the  acquisition  of  both  types  of 
data  at  the  same  x station  on  the  cone  surface. 

Calculations  were  performed  at  the  following  conditions: 

M - 7.95  Re  - 4.2  x 105  (x  - 4.0  in.) 

x 

P.  • 259.3  psia 

Coo 

T - 1360°R  T /T  - 0.41 

too  w t°° 

a - 0°,  8°,  12°,  20°,  20.5°,  24°  Pr  - 0.72 
The  finite  difference  mesh  contained  48  x 50  points  except  for  the 
a ■ 24°  case  in  %iiich  a 48  x /0  mesh  was  used.  The  circumferential 
(4>)  step  size  was  4°  for  all  runs  with  the  0 step  size  varying  depend- 
ing r?  the  angle  of  attack  and  number  of  mesh  points  used.  All  runs 
in  this  study  were  converged  to  the  fifth  significant  digit  of  the 
flow  quantities.  Computer  run  times  for  these  solutions  ranged  from 
1.4  hours  (a  ■ 12°)  to  1.6  hours  (a  • 24°)  on  the  CDC  CYBER  74 
(equivalent  in  CPU  speed  to  a CDC  6600) . The  number  of  time  steps 
required  ranged  from  1860  (a  ■ 8°)  to  2480  (a  ■ 12°) . The  fact  that 
the  highest  angle  of  attack  nm  did  not  take  the  largest  number  of 
time  steps  is  due  to  variance  In  the  actual  CFL  condition  at  which 
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each  calculation  was  made.  More  consistent  results  are  obtained 
through  use  of  the  stability  criteria  set  forth  in  Appendix  C.  In 
order  to  provide  comparison  of  both  flow  field  and  surface  quantities, 
comparison  with  experimental  data  Is  shown  for  surface  pressure  and 
pitot  pressure  surveys.  In  addition,  velocity  vector  plots  are  shown 
for  the  calculations  and  pictorial  displays  of  cross  flow  field  fea- 
tures are  shown  in  order  to  compare  the  overall  agreement  of  the 
calculation  with  experiment.  (These  data  were  also  reported  in  Ref  39). 

Figures  7 through  10  compare  the  numerical  results  of  the  calcu- 
lation with  the  surface  pressure  data  of  Tracy.  At  this  point  it 
should  be  noted  that,  unless  otherwise  stated,  the  discrete  points 
shown  for  Tracy's  data  were  obtained  from  continuous  curves  in  Ref  30 
through  use  of  a digitizer.  Figure  7 gives  this  comparison  for 
a ■ 8°.  The  calculated  value  can  be  seen  to  fall  just  at  the  edge  of 
the  symbols  for  the  entire  circumference  of  the  cone.  This  corresponds 
to  a local  difference  of  1.6Z  at  if  ■ 0°,  the  windward  centerline  of  the 
cone.  The  surface  pressure  comparison  for  a ■ 12°  is  given  in  Fig.  8. 
The  agreement  in  this  case  is  somewhat  better  near  the  leeward  center- 
line  (<j>  * 180°)  but  the  local  difference  has  increased  to  2.5Z  at  <j)  ■ 0. 
The  same  type  of  agreement  exists  for  the  a - 20°  (Fig.  9)  and  a ■ 24° 

(Fig.  10)  cases  with  the  local  difference  for  a - 24°  becoming  5Z  at 

<P  ■ 0°.  To  provide  a further  comparison,  the  surface  pressure  at  (f>  ■ 0° 

obtained  by  Lubard  and  Helliwell  (Ref  16)  for  the  a - 12°  case  is 

identical  to  that  obtained  lr  the  present  study.  The  first  possible 
reason  for  the  discrepancy  between  calculation  and  experiment  on  the 
windward  side  is  that  Tracy  noted  an  approximate  0.5°  uncertainty  in 
angle  of  attack  in  the  wind  tunnel  used  for  the  experiment.  In  order 
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to  test  this  possibility,  the  a ■ 20°  run  was  reconverged  at  a ■ 20.5°. 
The  surface  pressure  at  <j)  ■ 0°  regained  only  1/2  of  the  4.7Z  difference 
between  the  experiment  and  the  a ■ 20°  solution,  indicating  that  not 
ail  of  the  discrepancy  could  be  due  to  the  uncertainty  in  angle  of 
attack.  The  second  possibility  which  should  be  considered  is  the  fact 
that  the  surface  pressure  taps  used  for  the  experiment  were  somewhat 
large  in  relation  to  boundary  layer  thickness  on  the  windward  side. 

In  any  case,  an  error  band  of  at  least  5Z  is  present  in  most  experi- 
mental studies. 

Sample  circumferential  pitot  surveys  are  presented  in  Figs.  11 
through  13.  The  angle  in  the  6 direction  for  comparison  with  experi- 
ment is  defined  by 

8 - 9c  + tan-^j^rr]  (15) 

where  y is  the  normal  distance  above  the  cone  surface  in  inches.  The 
Reynolds  number  for  the  calculated  results  is  4.2  x 10^.  Since  the 
pitot  surveys  were  made  at  constant  height  above  the  cone  surface, 
they  will  in  most  cases  intersect  the  bow  shock  at  some  point  in  the 
circumferential  traverse.  This  appears  in  the  pitot  pressure  plots  as 

_3 

a sharp  rise  from  the  free  stream  value  (8.7  x 10  ) to  the  maximum 

value  for  each  trace.  Probe  effects  cause  the  character  of  the  rise 
of  experimental  pitot  pressure  to  appear  remarkably  like  that  of  the 
rise  of  pitot  pressure  through  the  captured  shock  waves  obtained  in 
the  present  study,  with  decrease  in  transition  slope  and  loss  of  peak 
pressure  in  the  calculations  caused  by  the  necessary  damping  present. 
Note  that  the  pitot  pressure  transition  centerpoint  (the  shock  location 
criteria  used  by  Tracy)  of  the  captured  shock  in  all  of  the  calculated 
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results  Is  virtually  the  same  as  would  be  determined  from  the  experi- 
mental pitot  survey  and  that  the  calculated  pitot  pressure  quickly 
returns  to  the  proper  magnitude  on  the  high  pressure  side  of  the  shock. 
For  the  a ■ 12°  case  (Fig.  11),  the  uppermost  survey  (y  ■ 0.4  in.)  is 
entirely  in  the  inviscid  flow  region  for  both  calculation  and  experi- 
ment. The  survey  nearest  the  cone  surface  (y  ■ 0.05  in.)  is  in  the 
viscous  layer  for  more  than  90°  of  the  circumference  of  the  cone. 
Agreement  with  experiment  for  both  of  these  surveys  is  good.  The 
calculated  center  survey  (y  - 0.2  in.)  departs  from  the  experimental 
value  near  the  lee  centerline  ($  - 180°) . This  departure  is  due  to 
the  smaller  viscous  layer  present  in  the  calculation  as  compared  to 
experiment . 

The  a - 20°  case  (Fig.  12)  has  essentially  the  same  degree  of 
agreement  between  the  calculated  results  and  experiment  with  the  excep- 
tion that  the  departure  near  the  lee  centerline  is  concentrated  in  a 
region  of  approximately  6°  (<J>  ■ 174°  - 180°).  Agreement  outside  this 
region  is  generally  good.  Figure  13  illustrates  that  essentially  the 
same  agreement  and  region  of  discrepancy  exist  for  the  a - 24°  case. 

The  pitot  pressure  discrepancy  near  the  lee  centerline  can  be 
evaluated  more  clearly  from  Figure  14,  which  compares  pitot  pressure 
along  the  lee  centerline  for  the  a ■ 24°  case.  The  discrepancy,  as 
noted  above,  extends  for  approximately  6°  to  either  side  of  the  lee 
centerline  for  this  angle  of  attack.  It  is  attributed  to  a combina- 
tion of  locally  large  pitot  pressure  drop  due  to  finite  model  nose 
radius  and  to  nose  effects  resulting  from  the  non-conical  region  at 
the  nose.  Any  persistent  experimental  nose  effects  would  tend  to  be 
concentrated  along  the  lee  centerline  as  they  are  swept  back  along  the 
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cone.  This  phenomenon  can  be  observed  in  the  experimental  study  of 

Stetson  (Ref  22)  through  examination  of  the  oil  flow  pictures  presented 

therein.  For  instance.  Fig.  12  of  Ref  22  shows  a near  conical  cross 

flow  separation  originating  at  approximately  16%  of  the  cones  length 

from  the  nose.  As  can  be  observed  from  the  oil  flow  traces  prior  to 

this  point,  the  flow  originating  at  or  near  the  nose  has  been  swept 

around  the  cone  to  the  lee  side  by  the  time  the  16%  point  has  been 

reached.  Since  th » flow  is  primarily  axial  on  the  lee  side  after  the 

cross  flow  separation  begins,  the  results  of  the  effective  blunting 

and  the  non-conical  region  near  the  nose  remain  at  or  near  the  lee 

centerline  as  the  flow  continues  along  the  cone. 

• • 

To  demonstrate  the  pressure  distribution  obtained  through  the 
shock  layer,  Fig.  15  shows  static  pressure  in  the  0 direction  at  60° 
intervals  around  the  cone  for  the  a - 12°  case.  The  elevation  above 
the  cone  surface  is  given  in  radians  in  this  and  subsequent  computer 
generated  plots.  As  noted  In  Appendix  8,  the  damping  was  not  tailored 
in  the  6 direction  for  these  runs,  so  excess  smearing  of  the  shock  is 
evident  on  the  free  stream  side.  No  appreciable  harm  is  seen  to 
result  from  this  smearing. 

The  velocity  vector  plots  in  Figs.  16  through  18  provide  an  indi- 
cation of  cross-flow  streamline  patterns  and  illustrate  the  relative 
magnitude  of  cross-flow  velocity  for  all  but  the  lowest  momentum 
regions.  The  reverse  flow  region  is  thin  in  the  a ■ 12°  case  (Fig.  16) 
with  the  vortical  singularity  (defined  for  this  study  as  the  cross 
flow  stagnation  point  toward  which  streamlines  converge)  occurring 
near  the  edge  of  the  viscous  layer  3.2°  above  the  cone  surface.  The 
cross-flow  separation  point  is  at  <J>  ■ 163°  with  the  experimental  point 
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occurring  at  $ ■ 156°.  This  difference  Is  attributable  to  lack  of 
resolution  of  the  reverse  flow  region  at  this  angle  of  attack  in  the 
present  study,  as  both  Lin  and  Rubin,  and  Lubard  and  Helliwell  obtained 
a separation  point  nearer  that  of  the  experimental  data.  This  was 
supported  in  the  present  study  by  the  fact  that  the  separation  points 
in  the  higher  a runs  (with  more  mesh  points  in  the  reverse  flow 
regions)  were  virtually  identical  to  those  of  experiment. 

At  a ■ 20°  (Fig.  17)  the  extent  of  the  reverse  flow  region  and 
the  position  of  the  cross-flow  separation  point  have  reached  limiting 
values,  as  no  substantial  changes  in  these  can  be  seen  in  the  a ■ 24° 
run  (Fig.  18) . This  was  noted  for  the  separation  point  in  both  experi- 
ment and  theory  by  Lin  and  Rubin  (Ref  15).  However,  the  extent  of  the 
viscous  region  and  the  location  of  the  vortical  singularity  (as  deter- 
mined in  the  present  study)  continue  to  change  as  a increases.  The 
vortical  singularity  is  near  the  edge  of  the  viscous  layer  for  all 
cases  in  the  present  study  in  which  it  is  lifted  off  the  cone  surface. 

Although  agreement  with  experiment  so  far  has  been  shown  to  be 
good,  perhaps  the  best  evaluation  of  the  validity  of  the  conical  flow 
assumption  for  engineering  solutions  can  be  made  through  a pictorial 
representation  of  the  flow  field  features.  Figures  19  and  20  compare 
the  overall  results  of  the  present  study  with  the  calculation  of  Lubard 
and  Helliwell  and  Tracy's  data.  As  shown  in  Fig.  19,  the  position  of 
the  bow  shock  and  viscous  layer  edge  of  the  present  study  for  this 
case  are  essentially  those  obtained  by  Lubard  and  Helliwell.  Both 
numerical  studies  result  in  a shock  position  too  near  the  cone  and  a 
thinner  viscous  layer  at  this  a when  compared  with  experiment.  This 
is  considered  to  be  primarily  due  to  the  failure  to  account  for  nose 
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tip  effects  as  noted  above.  A weak  supersonic  region  is  present  in 
the  cross  flow  plane  and  is  shown  by  the  solid  lines  between  viscous 
layer  and  bow  shock.  The  vortical  singularity  position  is  shown  by 
the  open  oval  symbol  on  the  lee  centerline.  The  structure  of  the 
a ■ 24°  case  is  more  interesting  (Fig.  20).  Supersonic  cross  flow  is 
seen  to  exist  over  most  of  the  field  and  is  terminated  by  a complex 
sonic  line/shock  wave.  The  position  of  this  line  is  nearer  the  lee 
centerline  than  shown  by  Tracy.  The  small  amplitude  oscillations 
present  on  the  Imbedded  sonic  line/shock  as  it  nears  the  bow  shock 
can  be  seen  (through  careful  examination  of  Fig.  18)  to  be  caused  by 
oscillations  (in  velocity)  propagating  through  the  supersonic  cross 
flow  region  from  the  vicinity  of  the  bow  shock.  The  <J>  grid  spacing 
used  in  the  present  study  did  not  give  sufficient  resolution  to  deter- 
mine the  existence  of  the  lambda  portion  of  the  imbedded  shock,  as 
shown  by  Tracy.  However,  its  signature  may  be  surmised  in  the  hump 
present  in  the  calculated  viscous  layer  edge  at  an  angle  of  approxi- 
mately 30°  from  the  lee  centerline.  The  curious  flat  top  present  in 
the  calculated  viscous  layer  edge  as  compared  to  experiment  is  again 
attributed  to  the  failure  of  the  numerics  to  model  non-conical  nose 
effects. 

In  summary,  the  present  technique  is  seen  to  model  all  features 
present  in  the  experimental  study,  with  the  exception  of  the  non- 
conical  nose  effects.  The  computer  times  are  not  excessive  since  they 
are  based  on  a very  stringent  convergence  criteria  and  represent  an 
increase  in  time  per  mesh  point  calculated  of  only  40%  over  a compar- 
able inviscld  computer  program. 
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Stetson's  Data 


Stetson  (Ref  22)  conducted  an  experimental  study  using  a 5.6° 
half  angle  sharp  (0.001  In.  nose  bluntness)  cone  in  M ■ 14.2  flow. 

The  wind  tunnel  working  fluid  was  air  with  total  pressure  of  1600 
psla  and  a total  temperature  of  2050°R.  The  resulting  free  stream  Re 
is  given  as  0.62  x 10^/Ft.  Data  were  taken  at  angles  of  attack  of  0° 
to  14°  in  2°  increments  (plus  a run  at  5°) . Measurements  of  surface 
pressure  were  made  at  a number  of  stations  along  the  cone,  and  pitot 
surveys  were  made  perpendicular  to  the  free  stream  velocity  at  x/L  - 
0.87  for  a - 10°. 


. calculation  was  made  at 

the  following  conditions: 

[ - 14.2 

Re  - 7.9  x 105 

X 

(x/L  -0.75, 

* 1600.0  psia 
too 

L - 15.37  in.) 

- 2050°R 

T /T  - 0.29 

too 

W too 

Pr  -0.72 

The  finite  difference  mesh  contained  48  x 50  points  with  the  circum- 
ferential step  size  being  4°  as  for  the  previous  calculations.  The 
computer  run  time  (with  convergence  criteria  as  before)  was  2.2  hrs. 
of  CDC  6600  time  for  2700  time  steps.  In  order  to  provide  comparison 
of  both  flow  field  and  surface  quantities,  plots  of  calculated  results 
versus  experimental  data  are  shown  for  surface  pressure  and  selected 
cuts  through  the  pitot  pressure  surveys.  Also,  a velocity  vector  plot 
of  the  calculated  results  is  shown.  A plot  of  surface  pressure 
compared  with  Stetson's  data  comprises  Fig.  21.  The  agreement  with 
experiment  is  quite  reasonable  with  the  maximum  difference  of  approx- 
imately 16%  (an  compared  with  local  static  pressure)  occurring  in  the 
vicinity  of  the  pressure  minimum.  Pressure  scale  expansion  in  Fig.  21 
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and  Che  low  absolute  magnitude  of  pressure  in  this  region  make  the 
agreement  appear  less  adequate  than  actually  exists.  The  complete 
circumferential  pressure  distribution  is  qualitatively  comparable  to 
that  of  Fig.  9.  The  surface  pressure  scale  of  Fig.  21  is  equivalent 
to  0.0  to  0.6  on  Fig.  6.  The  calculated  cross  flow  separation  point 
(26°  from  lee  centerline)  agrees  quite  well  with  the  experimental 
value  (27°  from  lee  centerline) . 

Pitot  pressure  comparisons  at  three  $ stations  on  the  lee  side 
of  the  cone  are  given  in  Figs.  22  through  24.  The  first  comparison  at 
<j>  * 180°  (Fig.  22)  shows  a large  difference  in  position  (approximately 
27Z  at  the  bow  shock)  between  given  levels  of  pitot  pressure  for  the 
calculation  and  experiment.  Although  the  same  phenomenon  was  noted  in 
Tracy's  a - 24°  case  (Fig.  14),  the  difference  was  small  at  the  edge 
of  the  viscous  layer  and  the  calculated  bow  shock  was  displaced  very 
little  from  the  experimental  position.  However  in  the  present  case, 
the  bow  shock  is  displaced  nearly  as  much  as  the  edge  of  the  viscous 
layer.  This  displacement  can  be  seen  to  persist  (although  much  smaller 
in  magnitude)  through  $ * 160°  (Fig.  23)  until  by  * 140°  (Fig.  24), 
the  agreement  is  nearer  that  expected  from  previous  results.  In  all 
three  plots  of  pitot  pressure,  the  trends  as  shown  by  the  calculation 
are  essentially  correct.  The  remaining  question  then  concerns  the 
possible  reason  for  the  discrepancy  between  the  calculated  and  experi- 
mental viscous  layer  thickness. 

Several  possibilities  for  the  cause  of  the  viscous  layer  discrep- 
ancy should  be  examined.  It  should  first  be  noted  that  the  perfect 
gas  assumption  results  in  a somewhat  higher  Reynolds  number  than  is 
quoted  in  Ref  22  (0.83  x 10**  vs.  0.62  x 10**).  Since  this  would  result 
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in  a somewhat  thinner  viscous  layer  for  the  calculation,  the  Reynolds 
number  was  lowered  20%  and  the  run  was  reconverged.  The  change  in 
shock  position  on  the  leeward  side  was  a maximum  of  0.25°  (as  might 
have  been  expected  for  this  small  change  in  Re) . The  resulting  small 
change  in  shock  position  indicated  that  a change  in  Re  of  at  least  an 
order  of  magnitude  would  be  necessary  for  better  agreement.  A more 
likely  candidate  for  at  least  part  of  the  discrepancy  in  viscous  thick- 
ness is  the  effective  viscous  nose  blunting  present  in  the  experiment. 
This  is  clearly  shown  in  Fig.  25  (taken  from  Ref  22)  and  would  account 
for  at  least  one-fourth  of  the  discrepancy  (assuming  no  change  in  bow 
shock  angle,  which  of  course  could  not  be  guaranteed). 

The  value  of  ;<  for  this  calculation  was  3.22.  Without  further 
information  concerning  the  source  of  the  viscous  layer  discrepancy, 
this  higher  value  of  x must  be  viewed  as  an  indication  that  the 
discrepancy  will  be  encountered.  (According  to  Fig.  6,  divergence  of 
the  present  technique  from  weak  interaction  theory  was  noted  at 
X ■ 1.55  for  Tracy's  conditions). 

To  illustrate  the  features  of  the  flow,  a plot  of  cross  flow 
velocity  vectors  comprises  Fig  26.  This  plot  illustrates  that  the 
lee  side  vortices  and  the  vortical  singularity  which  existed  in  Tracy's 
cases  exist  here  also.  The  vortical  singularity  occurs  near  the  edge 
of  the  viscous  layer  as  before. 

In  summary,  the  calculation  is  seen  to  model  the  features  and 
trends  present  in  Stetson's  experimental  study  with  the  only  apparent 
difficulty  being  the  failure  to  correctly  predict  lee  side  viscous 
layer  thickness  and  shock  position.  The  maximum  difference  of  approx- 
imately 2731  occurred  at  the  lee  centerline  (<(>  ■ 180°) . 
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McElderry  * s Data 


McElderry  (Ref  31)  analyzed  experimental  data  taken  for  a 6° 
half  angle  sharp  cone  38.06  In.  In  length  in  M ■ 6.05  flow.  It  should 
be  noted  that  these  data  also  appear  In  Ref  32  by  Rhudy  and  Baker. 

The  wind  tunnel  working  fluid  was  air  with  total  pressure  of  280.0  psla 
and  a total  temperature  of  850°R.  The  resulting  free  stream  Re  was 
given  as  5.0  x 10^/Ft.  Data  were  taken  at  angles  of  attack  of  0°  to 
12°  In  3°  Increments.  Surface  pressure  measurements  were  made  at 
several  locations  along  the  cone  with  pitot  surveys  made  perpendicular 
to  the  cone  centerline  at  x/L  -0.75  and  0.97.  Heat  transfer  data  were 
taken  on  a thin  walled  model  and  were  used  to  locate  the  boundary 
layer  transition  region  on  the  cone.  References  to  transition  contained 
In  the  remainder  of  the  discussion  of  these  data  are  based  on  the 
transition  location  analysis  presented  in  Ref  31.  At  the  experimental 
free  stream  Reynolds  number  of  5 x 10^/Ft. , transition  was  essentially 
complete  on  the  lee  side  by  x/L  ~ 0.3.  Since  the  surface  pressure 
measurement  nearest  the  tip  was  at  x/L  ■ 0.31,  nearly  all  lee  side 
pressure  data  were  taken  in  the  turbulent  flow  region. 

The  first  attempts  made  to  calculate  this  case  were  at  a Reynolds 
number  of  15.0  x 10^  which  corresponds  to  the  fully  turbulent  flow 
region  near  the  base  of  the  cone.  As  expected,  the  very  thin  boundary 
layer  present  at  this  Reynolds  number  created  extreme  difficulty  in  the 
vicinity  of  the  interaction  with  the  lee  side  imbedded  shock  wave 
(The  flow  at  this  Reynolds  number,  when  considered  without  provision 
for  turbulence  modeling,  could  be  characterized  as  nearly  inviscid  and 
therefore  contains  some  of  the  difficulties  encountered  with  the 
inviscid  solution) . No  combination  of  damping  and  step  size  was  found 
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which  would  lead  to  convergence  at  the  above  Reynolds  number.  The 
Reynolds  number  was  then  reduced  by  a factor  of  ten  to  1.5  x 10^  and 
the  run  was  successfully  completed.  This  reduced  Reynolds  number 
corresponds  to  a station  on  the  cone  just  prior  to  the  onset  of  transi- 
tion (as  determined  in  the  experimental  study)  at  the  noted  angle  of 
attack  and  therefore  was  in  a region  where  the  laminar  flow  equations 
were  still  valid. 

The  conditions  for  the  successful  run  were: 

M * 6.05  Re  ■ 1.5  x 10^  (x  s 4 in.) 

x 

- 280.0  psia  T /T  -(161 

too  w 

T - 850°R  Pr  - a 72 

too 

a - 12° 

The  finite  difference  mesh  contained  48  x 50  points  with  the  <j>  step 
size  being  4°  as  before.  Since  considerable  experimentation  with 
time  step  size  and  damping  factor  was  necessary  during  this  run,  the 
computer  time  of  2.3  hours  was  more  than  would  be  required  for  a rerun 
of  this  or  a similar  case. 

The  initial  behavior  of  this  solution  was  essentially  the  same  as 
that  for  the  higher  Reynolds  number  mentioned  above.  That  is,  the  lee 

side  Imbedded  shock  is  strong  until  very  near  the  cone  surface  (due  to 

the  thin  viscous  layer)  . This  resulted  in  a strong  compression  in  the 
vicinity  of  the  cone  surface  which  was  then  followed  by  an  over- 
expansion and  then  a recompression  to  the  cone  lee  centerline  point 
($  * 180°) . Although  this  appeared  to  be  a persistent  solution  to  the 
governing  equations  as  numerically  approximated,  the  overexpansion 
became  strong  enough  to  drop  pressure  (and  thereby  temperature)  below 
zero  which  then  destroyed  the  calculation.  Since  normal  stress  damping 
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was  found  to  be  effective  in  damping  starting  transients  in  the  lee 
side  region  for  Stetson's  data,  it  was  also  tried  here.  After  some 
experimentation,  it  was  found  that  the  overexpansion  could  be  controlled 
(and  eventually  made  to  disappear  entirely)  by  a damping  factor  of  225.0 
in  the  8 direction  and  100.0  in  the  (p  direction.  These  factors  are 
scaled  by  the  Reynolds  number  and  therefore  correspond  to  factors  of 
63.0  and  28.0  respectively  at  Tracy's  Reynolds  number.  In  contrast  to 
the  use  of  damping  to  control  starting  transients,  the  damping  in  this 
case  could  not  be  removed  as  the  overexpansion  would  then  reappear. 

The  above  damping  factors  were  then  maintained  until  the  solution  con- 
verged. Since  the  maintenance  of  normal  stress  damping  in  the  viscous 
layer  may  shift  the  density  distribution  through  the  layer  somewhat, 
changes  in  the  agreement  with  experiment  might  be  expected. 

Unfortunately,  experimental  measurements  were  not  made  far  enough 
forward  on  the  cone  to  provide  comparison  with  the  calculated  solution 
at  the  same  Reynolds  number.  The  first  surface  pressure  measurements 
were  at  x/L  - 0.31  at  which  point  transition  was  essentially  complete 
on  the  lee  side  at  12°  angle  of  attack.  The  dual  lee  side  separations 
associated  with  turbulent  flow  were  present  in  che  experiment  at  the 
x/L  *0.31  point  indicating  that  the  comparison  could  not  possibly  be 
exact  as  the  calculation  only  had  the  single  cross  flow  separation 
which  is  generally  associated  with  laminar  flow.  However,  this  cal- 
culation provides  a unique  opportunity  to  examine  the  adequacy  of  the 
conical  symmetry  assumption  for  turbulent  flow.  It  can  be  postulated 
that  if  the  conical  symmetry  assumption  is  to  be  valid  in  the  turbulent 
regime,  an  examination  of  the  flow  field  just  prior  to  the  start  of 
transition  would  differ  only  in  the  details  of  the  viscous  layer  when 
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compared  conically  with  a measurement  far  downstream  of  transition. 

To  test  this  hypothesis,  a comparison  of  the  calculation  Oc/L  * 0.10, 

Re  ■ 1.5  x 10^)  to  the  experimental  surface  pressure  (x/L  « 0.955, 

Re  ■ 15.0  x 10^,  Ref  31)  and  selected  pitot  surveys  (x/L  ■ 0.75, 

Re  ■ 12.0  x 10^,  Ref  32)  is  shown  in  Figs.  27  through  29.  In  this 

instance,  a conical  comparison  implies  examining  flow  quantities  at  two 
distinct  spherical  radii  (points  A and  B of  Fig.  1)  on  each  of  a series 
of  rays  passing  through  the  vertex  of  the  cone.  If  the  flow  quantities 
have  nearly  the  same  vector  and  scalar  values  at  points  A and  B on  a 

majority  of  the  rays,  then  the  flow  is  essentially  conical  in  charac- 

ter. In  the  present  comparison,  point  A corresponds  to  the  calculated 
values  at  a point  on  the  cone  upstream  of  transition  and  point  B corre- 
sponds to  the  experimental  measurements  downstream  rf  transition. 

The  experimental  surface  pressure  at  x/L  - 0.955  is  shown  conically 
compared  with  the  calculated  values  in  Fig.  27.  The  agreement  is  quite 
good  in  both  trend  and  absolute  value  until  the  vicinity  of  the  experi- 
mental primary  separation  at  <f>  ■ 126°  is  reached.  From  this  point  on 
toward  the  lee  centerline,  the  absolute  value  of  the  pressure  r.hows 
varying  agreement  with  experiment.  The,  calculated  separation  point 
(4>  ■ 163°)  is  near  the  experimental  secondary  separation  point 
(<t>  ■ 166°).  To  place  the  conical  comparison  in  perspective,  the  agree- 
ment of  the  present  calculation  with  experiment  is  far  better  than  is 
obtained  by  the  inviscld  techniques  shown  in  Ref  31. 

Pitot  pressure  surveys  for  both  calculation  and  experiment  normal 
to  the  cone  centerline  are  presented  in  Figs.  28  through  30.  The 
first  survey  is  at  ij)  ■ 180°  (Fig.  28).  Both  experiment  and  calculation 
are  characterized  by  very  high  speed  flow  prevailing  nearly  to  the  cone 
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surface.  For  example,  the  first  two  mesh  points  above  the  surface  In 
the  calculation  have  local  Mach  numbers  of  5.0  and  5.9  respectively. 

The  oscillations  occurring  between  1.0  and  2.0  inches  in  the  calcula- 
tion are  numerical.  However,  the  apparent  oscillation  in  the  calcula- 
tion occurring  between  0.0  and  1.0  in.  is  in  reality  the  solutions 
reflection  of  the  fact  that  a high  speed  "jet"  of  flow  occurs  (as 
evidenced  by  tba  high  local  Mach  numbers  noted  above)  on  the  lee  center- 
line  under  the  contra-rotating  vortices.  The  existence  of  this  jet  is 
confirmed  by  examination  cf  Fig.  31  which  slows  the  calculated  Mach 
number  distribution  along  lines  of  constant  <f>.  The  vor  :ex  structure 
is  seen  to  be  centered  at  $ * 170°.  The  local  Mach  number  clearly 
increases  at  the  first  two  mesh  points  above  the  surface  as  the  lee 
centerline  (<f>  - 180°)  is  approached.  The  apparent  oscillation  on  the 
<J>  - 180°  line  between  0.0  and  1.0  in.  is  in  reality  an  indication  of 
vortex  position.  Tils  region  of  high  speed  flow  can  also  be  readily 
seen  in  the  pitot  maps  in  Ref  32  at  x/L  ■ 0.97.  However,  the  resolu- 
tion of  the  lines  of  constant  pitot  pressure  is  insufficient  for  the 
x/L  ■ 0.75  surveys  in  Ref  32  to  be  able  to  plot  with  certainty  this 
"jet"  in  the  experimental  data  on  the  graph  shown  here.  The  pitot 
survey  at  <f>  ■ 160°  (Fig.  29)  shows  much  the  same  trend  as  the  <f>  ■ 180° 
survey  except  for  lower  momentum  flow  near  the  surface  and  a slightly 
thinner  viscous  layer  in  the  calculation  than  in  the  experiment.  At 
<j>  * 140°  (Fig.  30)  the  conical  agreement  of  the  calculation  with 
experiment  i.1  very  good. 

The  global  conical  agreement  with  experiment  can  be  examined 
through  the  characteristics  of  the  cross  flow  velocity  vector  plot  in 
Fig.  32.  One  Immediately  obvious  feature  of  this  plot  is  that  the 
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vertical  singularity  Is  not  lifted  off  the  surface  even  though  the  angle 
of  attack  is  twice  the  cone  hal.  angle.  This  is  in  relaity  evidence  of 
the  fact  that  the  viscous  layer  instead  of  having  a large  hump  on  the 
lee  side  now  has  a "cusp"  at  the  lee  centerline  with  a nearly  sharp  local 
minimum  in  viscous  layer  thickness.  The  experimentalist  have  observed 
this  behavior  of  the  viscous  layer  edge  and  it  has  been  noted  in  vapor 
screens  and  pitot  surveys  taken  in  turbulent  flow  on  cones  in  Refs  31, 

33,  and  34  (see  Fig.  33  taken  from  Ref  33).  Also,  comparing  the  extent 
of  the  lee  side  vortices  present  in  the  laminar  region  calculation  with 
the  pitot  maps  in  Refs  31  and  32,  it  is  apparent  that  the  essential 
size  and  character  of  the  primary  vortices  are  maintained  through 
transition  and  beyond. 

The  differences  in  the  characteristics  of  the  viscous  layer  and 
surface  pressure  distribution  with  increased  Reynolds  numbers  are 
graphically  and  plctorlally  shown  in  Fig.  33  which  is  taken  from  Ref  33. 
The  "cusp"  viscous  layer  structure  shown  in  Fig.  33  for  Rainblrd's 
experimental  results  has  been  considered  to  be  primarily  a turbulent 
flow  regime  phenomenon  on  conical  bodies.  However,  the  calculation 
made  at  McElderry's  conditions  in  the  present  study  demonstrates  that 
the  "cusp"  is  present  in  the  viscous  layer  upstream  of  transition.  The 
viscous  layer  structure  of  Rainbird  shown  in  Fig.  33  is  therefore  a 
phenomenon  occurring  in  the  laminar  flow  regime  which  persists  into  the 
turbulent  regime. 

In  summary,  the  solution  at  a location  on  the  cone  just  upstream 
of  the  boundary  layer  transition  location  is  shown  to  exhibit  features 
previously  associated  with  turbulent  flows  over  conical  bodies.  The 
agreement  of  this  solution  with  experiment  is  amazingly  good  when  it 
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is  compared  conically  with  experimental  measurements  made  far  down- 
stream in  the  turbulent  flow  regime.  The  adequacy  of  the  conical 
symmetry  assumption  for  turbulent  flow  is  therefore  assured  if  a 
reasonable  means  can  be  found  to  model  viscous  layer  turbulence. 


V.  Conclusions 


Inviscid  computations  for  hypersonic  cone  flow  have  proven 

to  be  unstable  at  angles  of  attack  greater  than  twice  the  cone 

half  angle.  In  addition,  previously  obtained  solutions  including 

the  effects  of  viscosity  have  required  input  of  experimentally 

obtained  data  and/or  have  not  been  demonstrated  for  a/6  >2.0. 

c 

In  order  to  resolve  these  difficulties,  solutions  have  been 
demonstrated  for  the  complete  flow  field  around  sharp  cones  in- 
cluding viscous  effects.  These  solutions  were  obtained  by  applying 
a conical  symmetry  assumption  to  the  complete  Navier-Stokes  equations 
and  numerically  integrating  the  resulting  equation  set  by  use  of 
MacCormack's  finite  difference  scheme.  Integrations  were  performed 
at  selected  sets  of  boundary  conditions  for  which  previously  pub- 
lished experimental  data  were  available.  The  solutions  obtained 
were  for  a Mach  number  range  of  6.05  to  14.2,  a local  Reynolds 
number  range  of  0.4  x 10^  to  1.5  x 10^,  and  angles  of  attack  from 
a/0c  **  0.0  to  2.4.  Stability  of  the  solutions  was  demonstrated 
through  use  of  a stringent  convergence  requirement.  A technique 
(normal  stress  damping)  was  developed  and  used  to  reduce  or  eliminate 
the  numerical  oscillations  which  developed  during  the  integration  in 
the  vicinity  of  strong  shock  waves.  The  results  of  the  integration 
were  compared  with  the  experimental  studies  of  Tracy,  Stetson,  and 
McElderry.  The  following  points  are  presented  concerning  the  conduct 
of  the  integration  and  the  comparison  with  experiment: 

(1)  Agreement  with  the  experimental  data  of  Tracy  was  gen- 
erally good.  The  major  features  of  the  experiment  (shock  waves, 
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sonic  lines,  and  viscous  layers)  were  present  and  essentially 
correct  in  location  and  magnitude  in  the  results  of  the  integra- 
tion. The  primary  discrepancies  noted  in  the  comparison  with 
Tracy's  data  were  an  approximate  5%  difference  in  surface  pres- 
sure near  the  windward  centerline  and  an  error  in  pitot  pres- 
sure in  a small  region  near  the  lee  centerline  for  the  ct/6^  ,>2.0 
cases.  The  latter  discrepancy  was  concluded  to  be  due  to  failure 
of  the  present  technique  to  account  for  non-conical  nose  effects. 
Agreement  of  the  results  of  the  integration  with  the  experimental 
data  of  Stetson  was  less  adequate  than  was  shown  to  exist  with  the 
data  of  Tracy.  The  agreement  with  Stetson's  surface  pressure 
measurements  was  quite  reasonable.  However,  the  position  of  the 
lee  side  shock  wave  in  the  results  was  27 1 nearer  the  cone  surface 
at  the  lee  centerline  than  in  the  experiment.  This  discrepancy  was 
concluded  to  be  primarily  caused  by  inability  of  the  present  tech- 
nique to  model  the  effective  nose  blunting  existing  in  the  experi- 
ment which  resulted  from  viscous  effects  at  this  high  Mach  number 
(14.2).  The  apparent  effect  of  this  nose  blunting  was  much  greater 
in  extent  in  Stetson’s  case  (hypersonic  similarity  parameter, 

X - 3.0)  than  in  Tracy's  cases  (x  a 1.0). 

(2)  The  present  technique  removes  the  angle  of  attack  limita- 
tions present  in  invlscid  calculations.  The  instabilities  associated 
with  the  point  where  the  high  strength  lee  side  imbedded  shock  wave 
reaches  the  surface  are  removed  by  the  Inclusion  of  viscous  effects 
in  the  physical  modeling  of  the  flow  field. 

(3)  The  adequacy  of  the  conical  symmetry  assumption  is  indi- 
cated for  the  turbulent  regime  on  conical  bodies.  A solution  obtained 
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just  prior  to  the  beginning  of  boundary  layer  transition  at  the  experi- 
mental conditions  of  McElderry  agrees  well  with  experiment  when  coni- 
cally projected  far  into  the  turbulent  regime. 

(4)  Normal  stress  damping  is  shown  to  provide  a physically 
based  means  to  provide  the  necessary  control  of  spurious  numerical 
oscillations  around  shock  waves  without  additional  computer  time.  It 
is  also  very  effective  for  control  of  overexpansion  and  for  control 
of  starting  transients  due  to  ill-suited  initial  conditions. 

In  summary,  the  present  technique  was  shown  to  be  a viable 
means  of  calculating  the  flow  over  cones  at  angle  of  attack  including 
viscous  effects.  Agreement  with  experiment  was  quite  good  for  x 1 
with  some  discrepancy  encountered  in  the  lee  side  viscous  layer  as  x 
Increased.  Angle  of  attack  limitations  encountered  with  the  inviscld 
equations  were  removed.  This  technique  should  therefore  be  considered 
as  an  alternative  to  the  use  of  the  inviscld  flow  equations  for  future 
calculations  of  flow  over  conical  bodies. 
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VI . Recommendations 


Three  primary  directions  are  seen  for  extension  of  the  work 
presented  In  this  report.  The  first  Is  the  confirmation  of  the  tech- 
nique for  more  general  conical  configurations,  such  as  delta  wings. 

The  second  is  the  extension  of  the  Reynolds  number  range  for  which  the 
calculations  can  be  accomplished  through  inclusion  of  a turbulence 
model  in  the  formulation  of  the  equation  set.  '*''e  third  is  the  inves- 
tigation of  the  adequacy  of  the  conical  symmetry  assumption  for  lower 
supersonic  Mach  numbers. 

A further  item  pointed  out  by  this  study  is  the  need  for  re- 
examination by  the  experimentalist  of  the  techniques  used  for  collection 
of  experimental  data.  As  the  numerical  computation  capability  improves, 
heretofore  undetected  discrepancies  in  data  collection  techniques  may 
be  pointed  out.  An  example  is  the  underestimation  of  surface  pressure 
on  the  windward  side  of  the  cone  in  the  present  study.  This  under- 
estimation in  past  lnvlscid  studies  was  attributed  to  failure  to  include 
viscous  effects.  However,  two  studies  (the  present  study  and  Ref  16) 
are  available  which  exhibit  this  underestimation  even  though  viscous 
effects  were  included  in  the  calculated  results.  A third  study  (Ref  14) 
matched  the  windward  surface  pressure  for  Tracy's  data  only  at  a much 
lower  Reynolds  number  than  existed  in  the  data.  This  is  considered  to 
be  a clear  indication  that  the  details  of  the  use  of  surface  pressure 
taps  where  thin  hypersonic  boundary  layers  exist  should  be  further 
studied.  This  also  points  out  the  need  for  cooperation  and  collabora- 
tion between  the  experimentalist  and  the  numericist  in  the  effort  to 
improve  overall  aerodynamic  prediction  techniques. 
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Fig.  1 Cross  Section  of  Super 
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Fig.  2 Coordinate  System  and  Computational  Mesh  Shown  in  Relation 
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Fig.  3 Cutaway  Drawing  Showing  Typical  Particle  Paths  for  Supersonic  Inviscid 
Flow  about  Cone  at  High  a 
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Fig.  4 Approximate  Crossflow  Streamline  Pattern  for  Cone 
at  High  a 
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of  Viscous  Layer  Edge  and  Spherical  Radius 


o 


5 

=3 

o 

_l 

< 

o 


< 

> 


>- 

cc 

o 


z O 
uj  oo 
00  — 
UJ  > 

a z 

Q_  — 


O | 


uo 

CM 


00 

£°IX  ’d 


VO 

DO- 

TH 

PH 


Comparison  of  Calculated  Cone  Surface  Pressure  with  that  Obtained  from  Hypersonic 
Weak  Interaction  Theory 
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Fig.  8 Comparison  of  Calculated  and  Experimental  Surface 
Pressure  at  Tracy’s  Conditions,  a - 12° 
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Fig.  9 Comparison  of  Calculated  and  Experimental  Surface 
Pressure  at  Tracy's  Conditions,  o - 20° 
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Fig.  11  Comparison  of  Circumferential  Pitot  Pressure  Surveys 
at  Constant  Height  Above  the  Cone  Surface  at  Tracy's 
Conditions,  a ■ 12° 
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Fig.  13  Comparison  of  Circumferential  Pitot  Pressure 
Surveys  at  Constant  Height  Above  the  Cone 
Surface  at  Tracy’s  Conditions,  a - 24° 
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Fig.  16  Cross  Flow  Velocity  Vectors  Plotted  as  Seen  on  Spherical 
Surface  at  Tracy’s  Conditions,  a - 12° 
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Fig.  17  Cross  Flow  Velocity  Vectors  Plotted  as  Seen  on 
Spl.erical  Surface  at  Tracy's  Conditions,  a ■ 20c 
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Fig.  18  Cross  Flow  Velocity  Vectors  Plotted  as  Seen  on 
Spherical  Surface  at  Tracy's  Conditions,  a - 24° 
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Fig.  19  Comparison  of  Calculated  and  Experimental 
Flow  Field  Features  at  Tracy's  Conditions, 
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Fig.  21  Comparison  of  Calculated  and  Experimental  Surface  Pressure 
at  Stetson's  Conditions 
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Fig.  22  Comparison  of  Pitot  Pressure  Survey  at  Constant  4> 
at  Stetson's  Conditions,  <J>  - 180° 


71 


DISTANCE  FROM  SURFACE  ( IN.  5 


Fig.  23 


Comparison  of  Pitot  Pressure  Survey  at  Constant 
<t>  at  Stetson's  Conditions,  $ * I60b 
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Fig.  24  Comparison  of  Pitot  Pressure  Survey  at  Constant 
<j>  at  Stetson’s  Conditions,  $ - 140® 
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Fig.  26  Cross  Flow  Velocity  Vectors  Plotted  as  Seen  on 
Spherical  Surface  at  Stetson's  Conditions 
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Fig.  27  Conical  Comparison  of  Calculated  (Rex  ■ 1.5 
x 10®)  and  Experimental  (Re^  ■ 15.0  x 10**) 
Surface  Pressure  at  McElderry's  Conditions 
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Fig.  28  Conical  Comparison  of  Calculated  (Re^j  ■ 1.5 
x 10°)  and  Experimental  Rex  * 15.0  x 10°) 
Pitot  Pressure  Survey  at  Constant  <p  at 
McElderry’s  Conditions,  $ ■ 180° 
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Fig.  29  Conical  Comparison  of  Calculated  (Re^j  “1*5 
z 10°)  and  Experimental  (Rejj  ■ 15.0  x 10°) 
Pitot  Pressure  Survey  at  Constant  <t>  at 
McElderry's  Conditions,  $ * 160° 
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Fig.  30  Conical  Comparison  of  Calculated  j[F.ex  ■ 1.5  x 10  ) 
and  Experimental  (Re^j  - 15.0  x 10°)  Pitot  Pressure 
Survey  at  Constant  $ at  McElderry's  Conditions, 
t - 140° 
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Fig.  31  Calculated  Mach  Number  Distributions  at 
Constant  <f>  at  McElderry's  Conditions 
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Fig.  32  Cross  Flow  Velocity  Vectors  Plotted  as  Seen  on 
Spherical  Surface  at  McElderry's  Conditions. 
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All  Measurements  Have  Been  Normalized  by  the  Corresponding 
Zero  Angle-of-Attack  Value. 


<K  = 180° 


8.0  ; a/0  =1.6 

H =0.81  ;R  =4x10* 


RAINBIRD 

M = 4.25  ; a/0  = 2.1 

OO 

H =0.03  ; x = 6x  10( 


Fig.  33  Comparison  of  Experimental  Laminar  and  Turbulent 
Viscous  Layer  Features  (Taken  from  Ref  33) 
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Fig.  34  Demonstration  of  the  Effect  of  Normal  Stress  Damping 
on  Calculated  Static  Pressure  Along  9,  a ■ 0° 
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APPENDIX  A 


Boundary  Conditions 

The  solution  of  the  Navier-Stokes  equations  by  numerical  methods 
is  in  reality  a boundary  value  problem  subject  to  the  same  require- 
ments and  constraints  as  the  solution  of  the  set  by  analytic  tech- 
niques. However,  the  mathematical  theory  concerning  existence  and 
uniqueness  of  solutions  for  this  equation  set  has  not  been  resolved 
(Ref  35).  This  statement  is  especially  true  for  compressible  flows. 
The  only  course  then  open  to  the  engineer  is  to  presume  a well  posed 
problem  exists  when  the  known  physically  derived  auxiliary  conditions 
are  imposed.  If  the  integration  of  the  equation  set  subject  to  these 
conditions  yields  a steady  state  solution  (which  of  course  is  not 
guaranteed) , then  experimental  results  for  these  same  conditions  must 
be  used  to  verify  that  the  solution  obtained  is  indeed  the  correct 
one.  When  this  is  done  for  several  related  sets  of  conditions  with 
good  results,  confidence  in  the  adequacy  of  the  solution  technique 
for  use  where  experimental  results  are  not  available  is  Increased. 
However,  the  user  of  the  technique  for  this  purpose  must  always  take 
care  to  examine  the  solutions  obtained  for  non-physical  features,  as 
it  is  never  possible  to  completely  define  initially  the  limits  of 
usefulness  of  an  equation  set  in  any  meaningful  sense.  The  above  pro- 
cedure was  followed  in  this  study.  The  rest  of  this  appendix  will 
present  the  auxiliary  conditions  used  and  the  manner  in  which  they 
were  Imposed . 

Examination  of  Fig.  2 reveals  that  the  finite  difference  mesh  is 
a spherical  surface  with  the  inner  edge  bounded  by  the  cone  and  the 
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outer  edge  free.  The  mesh  shown  is  360°  circumferentially  in  extent. 
Since  two  known  lines  of  symmetry  exist  at  <J>  ■ 0°  and  4>  » 180°,  the 
extent  of  the  surface  can  be  readily  reduced  by  one-half.  This 
reduction  can  be  made  without  increase  in  the  number  of  boundary 
values  necessary  to  the  calculation  by  specifying  an  additional  line  of 
mesh  points  in  the  9 direction  and  at  the  required  value  of  The 
results  of  each  sweep  during  the  integration  are  then  reflected 
symmetrically  to  these  lines  about  the  lines  of  symmetry.  This  causes 
the  integration  (at  the  line  of  symmetry)  to  perceive  that  a mirror 
image  of  the  solution  is  evolving  on  the  other  180°  of  the  mesh  (as  it 
would  be  if  the  sweep  were  actually  360°) . This  procedure  therefore 
requires  no  a priori  knowledge  of  the  flow  quantities  at  the  lines  of 
symmetry  and  implies  that  the  following  conditions  are  being  imposed 


at  these  lines  of 

symmetry: 

at  <(> 

- 0° 

and  <j>  ■ 

180° 

3u 

3$  “ 

0 

& ■ 

0 

Note  that  these  two  conditions 
are  imposed  only  at  one  line 

— 

3<J> 

0 

3e  . 

3<f> 

0 

of  symmetry.  They  are  then 
derivable  from  the  governing 
equations  at  the  other. 

w ■ 

0 

3T 

6<p 

0 

Two  boundaries  therefore  remain  which  must  be  dealt  with. 

Examination  of  the  six  equations  (the  Navier-Stokes  equations, 
the  continuity  and  energy  equations,  plus  the  equation  of  state)  as 
modified  reveals  that  all  three  velocity  components  and  temperature 
appear  as  second  derivatives  in  9 and  <j>  with  pressure  and  density 
appearing  as  first  derivatives  in  9 and  <J>.  A solution  should  then  be 
obtainable  through  the  specification  of  two  boundary  conditions  for 
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each  of  Che  velocity  components  and  temperature  and  one  each  for 
pressure  and  density.  (In  reality,  the  algebraic  equation  of  state 
eliminates  the  requirement  for  specifying  both  pressure  and  density.) 
This  was  done  in  the  present  study  in  the  following  manner: 
at  the  cone  surface  0 ■ 0 , $ ■ 0°-+]80° 
u • v ■ w ■ 0 


T - T (4>) 

W 


at  the  outer  boundary  (physically,  the  free  stream) 
0 ■ 0 outer,  0 ■ 0°-*-180° 


(40 

(40 


Woc  ($) 


a uniform  free 
stream  is  assumed 


Note  that  although  the  quantities  are  specified  along  a line  of 

constant  0,  the  4>  distribution  is  completely  specified  also.  It 

should  also  be  noted  that  the  total  interne1,  energy  was  used  in  place 

of  the  temperature  as  *a  dependent  variable.  The  connection  is  simply 
V2 

e ■ C T + -r—  . These  conditions  are  now  sufficient  to  allow  solution 
v 2 

of  the  problem. 

Since  pressure  and  density  at  the  cone  surface  remain  free  (and 
indeed  must  do  so  to  prevent  overspecification  of  the  boundary  condi- 
tions) , they  must  be  determined  during  the  integration  process.  This 
was  done  in  the  present  study  by  analytically  evaluating  the  0 
momentum  equation  at  the  cone  surface  conditions.  The  resulting 
equation  in  non-dlmensionalized  form  solved  for  -If-  is: 


1e  . i 

38  k sin6 


4f  cosQ 
3 1 Re 


+ sin0 


36  W J 36 


j8  sin8  3u  4 sin9  32v  2_  3 , 1 &w. 

3 Re  30  3 Re  302  " 3 30^Re  3<j>; 

+ i_A_  *±) 

3^Re  30' 

This  equation  is  approximated  by  one  aided  second  order  accurate 
finite  differences  about  the  surface  mesh  point  of  Interest.  The 
resulting  difference  equation  is  solved  for  the  surface  pressure,  the 
only  quantity  in  the  equation  which  has  not  been  updated  when  the 
equation  is  applied  after  a sweep  through  the  mesh  of  either  the 
predictor  or  the  corrector.  This  equation  can  be  represented  by: 


- 2A8 
3k  sln0 


(u,v,w) 


pj+i,n 


1 i 

3 Pj+2,n 


with  J » 1,  and  i and  n indicating  time  level  and  <j>  position  in  the 
mesh  respectively.  The  density  at  the  surface  is  found  via  the  equa- 
tion of  state  and  all  flow  quantities  are  then  known  or  have  been 
determined  at  the  cone  surface. 

To  complete  the  specification  of  the  auxiliary  conditions, 
initial  values  of  the  flow  quantities  are  required  at  all  points  in 
the  mesh  in  order  to  begin  the  calculation.  In  the  present  study, 
this  was  done  by  setting  the  mesh  to  free  stream  conditions  including 
the  mesh  point  above  the  cone  surface  with  true  surface  boundary  condi- 
tions plus  free  stream  pressure  and  density  imposed  at  the  cone  surface. 
The  calculation  was  then  begun  which  resulted  in  a so  called  "impulsive" 
start  up.  (i.e.,  the  cone  is  brought  from  rest  to  free  stream  velocity 
instantaneously) . The  outer  and  inner  boundary  conditions  were  main- 
tained throughout  the  calculation. 
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In  summary,  a means  of  reducing  the  field  of  calculation  by  half 
through  symmetry  has  been  discussed.  The  requirements  for  boundary 
conditions  were  shown  and  the  manner  in  which  these  requirements  were 
satisfied  in  the  present  study  was  presented.  The  use  of  the  8 
momentum  equation  to  obtain  the  cone  surface  pressure  was  demonstrated 
and  the  initial  conditions  utilized  were  also  discussed. 
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APPENDIX  B 


Normal  Stress  Damping 

Previous  studies  (Refs  13,  25)  using  shock  capturing  finite 
difference  techniques  have  encountered  oscillations  In  the  vicinity  of 
strong  shock  waves  which  can  cause  solution  instabilities.  To  overcome 
this  difficulty,  artificial  terms  have  been  added  to  the  governing 
equations  to  provide  necessary  damping.  In  the  present  study,  using 
the  viscous  equations  and  shock  capturing,  oscillations  were  again 
encountered.  This  indicated  that  the  natural  viscous  terms  are  inade- 
quate to  overcome  the  series  truncation  error  at  the  existing  mesh 
spacing.  To  damp  these  oscillations  without  additional  programming 
complexity,  the  normal  stress  terms  (which  are  in  general  large  only 
near  shock  waves)  were  altered  by  increasing  the  second  coefficient  of 
viscosity  (A>^  - 2/3y)  . This  resulted  in  elimination  or  reduction  of 
the  oscillations  where  desired.  In  addition  to  improving  the  shock 
structure,  it  was  found  that  normal  stress  damping  was  extremely 
effective  in  removing  instabilities  caused  by  starting  transients  which 
result  from  ill-suited  initial  conditions.  The  details  of  implementa- 
tion and  conditions  for  proper  use  of  this  technique  are  given  below. 

As  shown  in  Eq  (4),  the  normal  stress  terms  contain: 

°o  • XRe  <V-”  ■ - fRe[  2“  81,10  + ?F(v  8lne>  + !?]  CB1) 

Normal  stress  damping  results  when  2u  sln6  is  removed  and  the 

remainder  of  the  term  is  multiplied  by  an  empirically  determined  factor 

8,  where  6 is  negative  (i.e.  $a  ) . (A  positive  8 does  not  produce  the 

n 

damping  effect  and  makes  the  calculation  unstable  except  near  8 ■ 1)  • 
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The  damping  terms  occur  both  In  the  energy  equation  and  In  the  H" 
matrix  resulting  from  the  coordinate  transformation.  It  was  found  that 
the  most  accurate  solutions  resulted  from  use  of  the  damping  multiplier 
in  the  derivative  term  only  (F_  or  G.)  of  the  momentum  equation  most 

W $ 

nearly  normal  to  the  shock  wave  for  which  damping  is  desired.  Inclu- 
sion of  the  multiplier  in  the  i?"  matrix  or  the  energy  equation  produced 
some  displacement  of  the  solution.  Also,  it  should  be  noted  that  in 
the  steady  limit, 


a ■ f(V*V)  - f 
n p 


(B2) 


(through  use  of  the  continuity  equation)  . This  implies  that  normal 
stress  damping  cannot  be  used  In  the  vicinity  of  compressible  flow 
boundary  layers  without  regard  for  possible  changes  in  the  density 
gradients. 

An  example  of  the  use  of  normal  stress  damping  is  given  in  Fig. 

34  for  axisymmetric  cone  flow  as  calculated  by  the  technique  of  this 
report.  The  plot  is  a comparison  of  the  calculated  static  pressure  in 
the  0 direction  at  Tracy's  free  stream  conditions  with  and  without  the 
inclusion  of  normal  stress  damping.  In  this  case  the  damping  is 
tailored  in  the  9 direction  by: 

B - 1 - (fe>  A1  - <fc>  *2  (B3) 


where  the  pressure  ratio  multipliers  were  empirically  determined  to  be 
A1  - 9 and  A2  ■ 2.  The  use  of  the  damping  multiplier  (Fig.  34)  resulted 
in  the  reduction  of  the  numerical  oscillations  near  the  shock  wave. 
Although  a small  offset  in  pressure  occurs  outside  the  boundary 
laye",  the  surface  pressure  and  shock  position  do  not  change 
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appreciably  when  damping  is  applied. 

For  Tracy '8  cases,  the  damping  multiplier  6 was  unity  over  the 
lee  40°  (<£  • 140°  - 180°)  of  the  flow,  where  viscous  effects  predomi- 
nate. It  was  then  varied  smoothly  around  the  cone  from  <p  ■ 140°  in 
accord  with  the  equation: 


8 . - 1 - 8' 

<*> 


p 

s 

A P 

<f>  - s 

<t>  ■ 140° 

p 

s 

-e- 

i 

o 

o 

- PsU  - 140° 

(B4) 


where  s indicates  body  surface  values  and  8"  was  empirically  varied 
from  15  (a  ■ 8°)  to  130  (a  ■ 24°).  8 was  constant  in  the  0 direction. 

For  Stetson's  case,  the  very  strong  lee  bide  expansion  during 
start  up  made  it  necessary  to  use  normal  stress  damping  there. 

However,  the  damping  in  this  region  was  removed  after  the  solution 
development  was  sufficient  for  the  expansion  to  become  less  strong. 

For  the  rest  of  the  field  (<j)  ■ 0°-KL40°)  , a more  sophisticated  tailoring 


of  the  damping  multiplier  was  used: 


col 


6 9,0  " 1 ' 


Psl»  ~PsU  - 140° 

[Psl((i  - 0°  - Ps  | <f>  ■ 140° 


l<t>  ■ oc 


PU  - 140°  J 


co2 


(B5) 


The  value  of  col  was  1.6  and  of  co2  was  1.1.  This  equation  was  based 
on  experiments  with  the  use  of  normal  stress  damping  which  showed  that 
the  effects  could  be  improved  by  tailoring  according  to  the  local 
value  of  static  pressure.  An  Improvement  was  noted  in  the  damping 
effect,  but  it  was  found  that  this  equation  did  not  always  increase 
the  magnitude  of  8 quickly  enough  during  initial  start  up.  This  was 
overcome  by  using  a high  initial  magnitude  of  8 for  approximately  the 
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first  100  tine  steps  of  the  run.  The  above  equation  was  used  for 
McElderry's  case  also.  However,  the  value  of  3 necessary  on  the  lee 
side  was  more  negative  than  the  output  of  the  above  equation.  The  high 
initial  magnitude  was  therefore  used  everywhere  in  the  flow  field. 

In  order  that  others  who  use  normal  stress  damping  may  benefit 
from  experience  gained  in  the  present  study,  the  following  observations 
and  rules  concerning  this  technique  are  offered: 

(1)  The  most  valuable  and  recommended  use  of  normal  stress  damping 
is  to  reduce  starting  transients  from  ill  suited  initial  conditions.  No 
permanent  history  of  the  use  of  damping  for  this  purpose  is  seen  in  the 
flow  after  the  removal  of  the  damping  and  subsequent  convergence. 

(2)  Normal  stress  damping  has  been  used  to  allow  capturing  of 
shock  waves  with  p2/pi  as  high  as  32.  The  only  practical  limit  to 
shock  strength  which  can  be  captured  through  use  of  this  technique 
appears  to  be  number  of  mesh  points  available  for  the  smearing  to 
take  place  over. 

(3)  Tailoring  of  the  damping  based  on  local  static  pressure  was 
found  to  give  reasonable  results.  Tailoring  based  on  other  quantities 
(such  as  temperature)  gave  unnecessary  offset  in  the  viscous  layer. 

(4)  Large  shifts  of  surface  pressure  or  of  static  pressure  behind 
shocks  when  damping  is  applied  usually  indicate  that  the  magnitude  of 

@ is  too  high  in  that  locale. 

(5)  The  best  results  are  obtained  behind  the  shock  wave  when 
some  oscillations  are  allowed  to  remain  in  the  free  stream. 

(6)  Care  must  be  taken  to  insure  that  the  viscous  stability  limit 
is  satisfied  with  3 included  in  the  manner  shown  in  Appendix  C. 

(7)  The  items  noted  below  are  observed  to  cause  undue  displacement 
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of  the  solution  at  high  values  of 

(a)  Inclusion  of  the  B multiplier  In  terms  resulting  from 
coordinate  transformations  (l.e.,  the  matrix). 

(b)  Inclusion  of  $ In  the  energy  equation. 

The  above  observations  and  rules  should  not  be  considered  as  axiomatic 
but  are  offered  as  guides  for  the  use  of  this  technique  in  other 
compresssible  flow  problems. 

Normal  stress  damping  is  a physically  based  technique  which  has 
proven  beneficial  in  damping  out  starting  transients  and  smoothing 
numerical  oscillations  in  the  vicinity  of  shock  waves.  The  3 contouring 
equations  and  levels  used  for  the  various  cases  were  presented.  Rules 
and  observations  were  also  given  to  guide  future  use  of  this  technique. 


APPENDIX  C 


Stability  Criteria 

In  order  to  maximize  computational  efficiency,  criteria  should 
be  used  to  determine  the  maximum  permissible  time  step  (At)  as  the 
integration  progresses.  As  noted  in  Section  III,  criteria  were  not 
used  during  the  calculations  presented  elsewhere  in  this  report. 

Instead,  At  was  fixed  during  each  complete  integration  at  Tracy's  condi- 
tions and  was  fixed  for  the  later  parts  of  the  integrations  at  Stetson's 
and  McElderry's  conditions.  The  value  at  which  At  was  fixed  was  at 
first  determined  by  increasing  At  until  instability  was  encountered 
during  the  starting  transient  phase  of  the  run.  The  run  could  then  be 
continued  to  convergence  by  reducing  At  10Z  to  15Z  from  the  unstable 
value.  As  experience  was  gained  with  the  technique  (and  more  was 
learned  about  the  effect  of  normal  stress  damping  on  stability),  it  was 
found  that  for  Tracy's  conditions,  stable  values  of  At  could  be  set 
initially  without  resorting  to  trial  and  error.  These  values  were  in 
general  O.6A0  to  O.7A0.  The  above  procedure  did  not,  however,  allow 
for  any  possible  increase  in  allowable  At  after  the  large  transient  phase 
of  the  calculation.  Difficulty  in  determining  allowable  At  was  also 
encountered  for  Stetsorfs  and  McElderry's  conditions  due  to  use  of  normal 
stress  damping  in  predominately  viscous  regions.  To  correct  these  two 
deficlences,  stability  criteria  for  use  during  the  integration  were 
postulated  and  confirmed  by  a repeat  of  previous  calculations.  The 
remainder  of  this  appendix  will  present  this  procedure  and  the  observed 
changes  in  computational  efficiency. 
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A complete  stability  analysis  for  MacCormack's  method  as  applied 
to  the  Navier-Stokes  equations  has  not  been  accomplished.  However,  use- 
ful estimates  of  allowable  step  sizes  can  be  made  through  the  use  of 
linear  theory  and  physical  arguments.  Claiming  that  the  physical  propa- 
gation of  information  (i.e.  fluid  signals)  cannot  outstrip  the  numerical 
propagation  results  in  the  following  Courant-Friedrich-Lewy  (CFL) 
condition  (Ref  36) : 

..  . 1.0 

At  £ r 

C I vj_  . I w|  CA92  + (sinB  A4>)2]  * 

A0  sine  A$  A0(sin6  At)  (Cl) 

This  expression  essentially  determines  allowable  At  in  regions  where 

viscous  effects  are  negligible.  The  remaining  task  is  to  determine  the 

effects  on  stability  of  the  stress  terms  and  of  normal  stress  damping. 

The  maximum  stable  At  in  viscous  regions  is  normally  a function  of 
the  smallest  physical  step  size  there.  In  the  present  study  A0«Acf> 
for  all  cases  calculated.  This  implied  that  the  stress  terms  in  the  0 
coordinate  direction  were  the  main  determinant  of  allowable  At  and  there- 
fore  that  the  stability  analysis  due  to  the  stress  terms  was  essentially 
a one  dimensional  problem.  In  the  case  of  finite  difference  operator 
time  splitting,  the  stability  analysis  also  has  been  shown  to  be  a one 
dimensional  problem  (Ref  25) . An  analysis  similar  to  that  of  Ref  25 
was  therefore  accomplished  and  is  presented  here. 

Equation  (2)  is  first  written  in  non-conservation  form: 

[L][Z]  t + [«J[Z]  e + [N][Z]  t + [«][Z]  e0  + [Q][z]  # 

» 

+ 0}(X)e<,  + [xirzj  + [«']■»  <“> 
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where  [ Z ] 


u and  the  subscripts  indicate  differentiation. 


The  remaining  matrices  will  be  defined  as  they  are  used.  The 
final  equation  results  when  [l]  ^ is  found  and  the  coefficient  matrices 

[L  ] ‘ [mJ  , etc.  are  formed  By  declaring  the  coefficient  matrices 
to  be  locally  constant,  it  can  be  shown  that  (Ref  37,  38)  when  each 
dimension  is  considered  separately  V £ 1.0  (lnvlscid  effect)  and 
v (viscous  effect) . The  symbol  v is  in  each  case  the  maximum 

eigenvalue  of  the  appropriate  coefficient  matrix  (for  example  [ij 
[m]  ).  Equation  (Cl)  above  accounts  for  the  lnvlscid  effect  on  stabil- 
ity so  that  it  is  only  necessary  to  consider  the  second  inequality 
and  obtain  values  of  v for  the  coefficient  matrices  of  the  second  deriv- 
ative terms  involving  0 in  Eq  (C2)  . 

The  matrix  f l]  ^ [W]  is  as  follows: 

[L] -1  [W  ] -f0  0 0 0 0 “1 

0 ^ 0 ex  0 0 

0 0 — 0 0 

pRe 


2p2RePr 


pRe  Pr 


The  maximum  eigenvalues  of  this  matrix  are  y/(pRe  Pr)  and  (2  + — ) /(pRe), 
The  second  eigenvalue  is  Included  due  to  the  undetermined  (and  sometimes 
large)  value  of  the  damping  coefficient  8. 
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Arguments  must  now  be  made  regarding  contribution  of  the  0<J>  and 


cross  derivative  terms  in  the  equation  set.  The  gofl  is  to  demonstrate 
that  the  elgnevalues  of [i/]  L[w]give  a minimum  At,  so  that  the  6<t>  and  <j>6 
terms  can  be  dropped  from  consideration.  First  assume  that  the  contri- 
bution of  these  terms  is  v £ y . This  effectively  increases  the 
contribution  (results  in  a lover  At)  of  these  terms  since  A0«A $.  As 
an  example,  the  matrix [L]”1[x]is 

Ortx]  - | 0 0 0 0 


0 

0 

0 
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0 

0 

0 
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0 
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pRe 

0 


&/2 

Re 
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0 

0 

0 
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(C4) 


The  maximum  eigenvalue  of  this  matrix  is  (0X/y)  V(pRe)  and  is  clearly 
less  than  the  eigenvalue  due  to  normal  stress  damping  of  the[L]  1[w] 
matrix.  Since  the[Lj  ^[sjmatrix  is  similar,  it  can  be  concluded  that 
sufficient  information  is  known  to  postulate  a complete  stability  cri- 
teria. By  substituting  the  eigenvalues  into  the  appropriate  stability 
criteria  and  labeling  each  resulting  At  as  to  its  source,  the  following 
complete  stability  criteria  is  postulated: 


At  * m min  (At  , At  , , At  . ) 

cfl  nsd’  vise 


(C5) 


where  m _<  1.0,  Atc^  ■ min  of  results  from  Eq.  (Cl) 

A . . pRe  A92 

4tnsd  ' 1 2 | (2+BX)| 

At  . - min  P?‘.  .P-r  M2 

vise  2 y 

where  min  implies  the  minimum  value  of  the  quantity  found  by  searching 
all  mesh  points  in  the  flow  field.  The  factor  m is  theoretically 
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limited  to  a value  of  1.0.  However,  uncertainties  and  approximations 
In  the  stability  analysis  sometimes  result  in  an  actual  maximum  value 
of  m greater  or  less  than  1.0. 

The  above  procedure  was  tested  at  o - 0°  for  Tracy's  conditions 
with  good  results.  However,  the  real  Issue  was  whether  the  criteria 
would  be  adequate  for  high  a with  the  use  of  large  amounts  of  normal 
stress  damping.  Therefore  the  a - 20°  case  at  Tracy's  conditions  was 
repeated  with  the  above  criteria  for  determining  At  in  continuous  use. 
Use  of  the  criteria  reduced  the  number  of  time  steps  necessary  for 
convergence  to  the  fifth  significant  digit  from  1700  to  960.  The 
factor  m for  this  run  was  1.0  for  the  first  300  time  steps  and  1.2 
tlereafter.  Instability  of  the  solution  occurred  for  m - 1.4.  The 
values  of  the  normal  stress  damping  coefficient  (6)  were  comparable 
with  those  used  in  the  original  run.  In  order  to  insure  that  the 
viscous  stability  limit  would  not  be  violated  when  using  very  high 
values  of  normal  stress  damping,  a run  was  made  at  Tracy's  conditions 
and  a - 0°  with  8 ■ 4000.  This  set  of  conditions  was  unstable  at  a 
value  of  At  determined  from  Eq  (Cl).  However  Vie  use  of  the  above 
criteria  reduced  the  value  of  At  sufficiently  to  allow  the  run  to  be 
continued  (although  at  a very  low  At) . 

In  summary,  stability  criteria  were  postulated  and  confirmed  which 
reduced  the  number  of  time  steps  required  for  convergence  by  44%  in  an 
example  calculation.  The  procedure  was  also  shown  to  account  for  the 
effects  of  normal  stress  damping  on  stability. 
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